38 #ifndef GETFEM_FOURTH_ORDER_H__ 
   39 #define GETFEM_FOURTH_ORDER_H__ 
   54   template<
typename MAT, 
typename VECT>
 
   57    const mesh_fem &mf_data, 
const VECT &A,
 
   61        "M(#1,#1)+=sym(comp(Hess(#1).Hess(#1).Base(#2))(:,i,i,:,j,j,k).a(k))");
 
   66     assem.
push_mat(
const_cast<MAT &
>(M));
 
   70   template<
typename MAT, 
typename VECT>
 
   71   void asm_stiffness_matrix_for_homogeneous_bilaplacian
 
   72   (
const MAT &M, 
const mesh_im &mim, 
const mesh_fem &mf,
 
   74     generic_assembly assem
 
   76        "M(#1,#1)+=sym(comp(Hess(#1).Hess(#1))(:,i,i,:,j,j).a(1))");
 
   80     assem.push_mat(
const_cast<MAT &
>(M));
 
   85   template<
typename MAT, 
typename VECT>
 
   86   void asm_stiffness_matrix_for_bilaplacian_KL
 
   87   (
const MAT &M, 
const mesh_im &mim, 
const mesh_fem &mf,
 
   88    const mesh_fem &mf_data, 
const VECT &D_, 
const VECT &nu_,
 
   90     generic_assembly assem
 
   91       (
"d=data$1(#2); n=data$2(#2);" 
   92        "t=comp(Hess(#1).Hess(#1).Base(#2).Base(#2));" 
   93        "M(#1,#1)+=sym(t(:,i,j,:,i,j,k,l).d(k)-t(:,i,j,:,i,j,k,l).d(k).n(l)" 
   94        "+t(:,i,i,:,j,j,k,l).d(k).n(l))");
 
   97     assem.push_mf(mf_data);
 
  100     assem.push_mat(
const_cast<MAT &
>(M));
 
  104   template<
typename MAT, 
typename VECT>
 
  105   void asm_stiffness_matrix_for_homogeneous_bilaplacian_KL
 
  106   (
const MAT &M, 
const mesh_im &mim, 
const mesh_fem &mf,
 
  107    const VECT &D_, 
const VECT &nu_,
 
  109     generic_assembly assem
 
  110       (
"d=data$1(1); n=data$2(1);" 
  111        "t=comp(Hess(#1).Hess(#1));" 
  112        "M(#1,#1)+=sym(t(:,i,j,:,i,j).d(1)-t(:,i,j,:,i,j).d(1).n(1)" 
  113        "+t(:,i,i,:,j,j).d(1).n(1))");
 
  117     assem.push_data(nu_);
 
  118     assem.push_mat(
const_cast<MAT &
>(M));
 
  135   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  146   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  147    const std::string &dataname1, 
const std::string &dataname2,
 
  159   template<
typename VECT1, 
typename VECT2>
 
  163     GMM_ASSERT1(mf_data.
get_qdim() == 1,
 
  164                 "invalid data mesh fem (Qdim=1 required)");
 
  190       s = 
"Grad_Test_u.(A*Normal)";
 
  192       s = 
"Grad_Test_u.(((Reshape(A,meshdim,meshdim)*Normal).Normal)*Normal)";
 
  194       s = 
"((Grad_Test_u')*A).Normal";
 
  197       s = 
"((((Grad_Test_u').Reshape(A,qdim(u),meshdim,meshdim)).Normal).Normal).Normal";
 
  199       GMM_ASSERT1(
false, 
"invalid rhs vector");
 
  200     asm_real_or_complex_1_param_vec(B, mim, mf, &mf_data, F, rg, s);
 
  203   template<
typename VECT1, 
typename VECT2>
 
  204   void asm_homogeneous_normal_derivative_source_term
 
  205   (VECT1 &B, 
const mesh_im &mim, 
const mesh_fem &mf,
 
  206    const VECT2 &F, 
const mesh_region &rg) {
 
  231     if (mf.get_qdim() == 1 && Q == 1)
 
  232       s = 
"Test_Grad_u.(A*Normal)";
 
  233     else if (mf.get_qdim() == 1 && Q == gmm::sqr(mf.linked_mesh().dim()))
 
  234       s = 
"Test_Grad_u.(((Reshape(A,meshdim,meshdim)*Normal).Normal)*Normal)";
 
  235     else if (mf.get_qdim() > 
size_type(1) && Q == mf.get_qdim())
 
  236       s = 
"((Test_Grad_u')*A).Normal";
 
  238              Q == 
size_type(mf.get_qdim()*gmm::sqr(mf.linked_mesh().dim())))
 
  239       s = 
"((((Test_Grad_u').Reshape(A,qdim(u),meshdim,meshdim)).Normal).Normal).Normal";
 
  241       GMM_ASSERT1(
false, 
"invalid rhs vector");
 
  242     asm_real_or_complex_1_param_vec(B, mim, mf, 0, F, rg, s);
 
  259   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  260    const std::string &dataname, 
size_type region);
 
  271   template<
typename VECT1, 
typename VECT2>
 
  272   void asm_neumann_KL_term
 
  273   (VECT1 &B, 
const mesh_im &mim, 
const mesh_fem &mf, 
const mesh_fem &mf_data,
 
  274    const VECT2 &M, 
const VECT2 &divM, 
const mesh_region &rg) {
 
  275     GMM_ASSERT1(mf_data.get_qdim() == 1,
 
  276                 "invalid data mesh fem (Qdim=1 required)");
 
  278     generic_assembly assem
 
  279       (
"MM=data$1(mdim(#1),mdim(#1),#2);" 
  280        "divM=data$2(mdim(#1),#2);" 
  281        "V(#1)+=comp(Base(#1).Normal().Base(#2))(:,i,j).divM(i,j);" 
  282        "V(#1)+=comp(Grad(#1).Normal().Base(#2))(:,i,j,k).MM(i,j,k)*(-1);" 
  283        "V(#1)+=comp(Grad(#1).Normal().Normal().Normal().Base(#2))(:,i,i,j,k,l).MM(j,k,l);");
 
  287     assem.push_mf(mf_data);
 
  289     assem.push_data(divM);
 
  294   template<
typename VECT1, 
typename VECT2>
 
  295   void asm_neumann_KL_homogeneous_term
 
  296   (VECT1 &B, 
const mesh_im &mim, 
const mesh_fem &mf,
 
  297    const VECT2 &M, 
const VECT2 &divM, 
const mesh_region &rg) {
 
  299     generic_assembly assem
 
  300       (
"MM=data$1(mdim(#1),mdim(#1));" 
  301        "divM=data$2(mdim(#1));" 
  302        "V(#1)+=comp(Base(#1).Normal())(:,i).divM(i);" 
  303        "V(#1)+=comp(Grad(#1).Normal())(:,i,j).MM(i,j)*(-1);" 
  304        "V(#1)+=comp(Grad(#1).Normal().Normal().Normal())(:,i,i,j,k).MM(j,k);");
 
  309     assem.push_data(divM);
 
  325   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  326    const std::string &dataname1, 
const std::string &dataname2,
 
  350   template<
typename MAT, 
typename VECT1, 
typename VECT2>
 
  354    const VECT2 &r_data, 
const mesh_region &rg, 
bool R_must_be_derivated, 
 
  356     typedef typename gmm::linalg_traits<VECT1>::value_type value_type;
 
  357     typedef typename gmm::number_traits<value_type>::magnitude_type magn_type;
 
  361     if (version & ASMDIR_BUILDH) {
 
  364         s = 
"M(#1,#2)+=comp(Base(#1).Grad(#2).Normal())(:,:,i,i)";
 
  366         s = 
"M(#1,#2)+=comp(vBase(#1).vGrad(#2).Normal())(:,i,:,i,j,j);";
 
  374       gmm::clean(H, gmm::default_tol(magn_type())
 
  375                  * gmm::mat_maxnorm(H) * magn_type(1000));
 
  377     if (version & ASMDIR_BUILDR) {
 
  379                 "invalid data mesh fem (Qdim=1 required)");
 
  380       if (!R_must_be_derivated) {
 
  383         asm_real_or_complex_1_param_vec(R, mim, mf_mult, &mf_r, r_data, rg,
 
  384                                         "(Grad_A.Normal)*Test_u");
 
  410   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  411    const std::string &multname, 
size_type region,
 
  412    const std::string &dataname = std::string(),
 
  413    bool R_must_be_derivated = 
false);
 
  430   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  431    const mesh_fem &mf_mult, 
size_type region,
 
  432    const std::string &dataname = std::string(),
 
  433    bool R_must_be_derivated = 
false);
 
  450   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  452    const std::string &dataname = std::string(),
 
  453    bool R_must_be_derivated = 
false);
 
  473   (model &md, 
const mesh_im &mim, 
const std::string &varname,
 
  474    scalar_type penalisation_coeff, 
size_type region, 
 
  475    const std::string &dataname = std::string(),
 
  476    bool R_must_be_derivated = 
false);
 
Generic assembly of vectors, matrices.
 
void push_data(const VEC &d)
Add a new data (dense array)
 
void push_mat(const MAT &m)
Add a new output matrix (fake const version..)
 
void assembly(const mesh_region ®ion=mesh_region::all_convexes())
do the assembly on the specified region (boundary or set of convexes)
 
void push_mi(const mesh_im &im_)
Add a new mesh_im.
 
void push_mf(const mesh_fem &mf_)
Add a new mesh_fem.
 
Describe a finite element method linked to a mesh.
 
virtual dim_type get_qdim() const
Return the Q dimension.
 
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
 
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
 
Describe an integration method linked to a mesh.
 
const mesh & linked_mesh() const
Give a reference to the linked mesh of type mesh.
 
structure used to hold a set of convexes and/or convex faces.
 
const mesh_region & from_mesh(const mesh &m) const
For regions which have been built with just a number 'id', from_mesh(m) sets the current region to 'm...
 
static mesh_region all_convexes()
provide a default value for the mesh_region parameters of assembly procedures etc.
 
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....
 
Model representation in Getfem.
 
void asm_stiffness_matrix_for_bilaplacian(const MAT &M, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT &A, const mesh_region &rg=mesh_region::all_convexes())
assembly of .
 
void asm_normal_derivative_dirichlet_constraints(MAT &H, VECT1 &R, const mesh_im &mim, const mesh_fem &mf_u, const mesh_fem &mf_mult, const mesh_fem &mf_r, const VECT2 &r_data, const mesh_region &rg, bool R_must_be_derivated, int version)
Assembly of normal derivative Dirichlet constraints  in a weak form.
 
void asm_normal_derivative_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg)
assembly of .
 
void asm_normal_source_term(VECT1 &B, const mesh_im &mim, const mesh_fem &mf, const mesh_fem &mf_data, const VECT2 &F, const mesh_region &rg)
Normal source term (for boundary (Neumann) condition).
 
size_t size_type
used as the common size type in the library
 
GEneric Tool for Finite Element Methods.
 
size_type add_bilaplacian_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname, size_type region=size_type(-1))
Adds a bilaplacian brick on the variable varname and on the mesh region region.
 
size_type add_bilaplacian_brick_KL(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname1, const std::string &dataname2, size_type region=size_type(-1))
Adds a bilaplacian brick on the variable varname and on the mesh region region.
 
size_type add_normal_derivative_Dirichlet_condition_with_multipliers(model &md, const mesh_im &mim, const std::string &varname, const std::string &multname, size_type region, const std::string &dataname=std::string(), bool R_must_be_derivated=false)
Adds a Dirichlet condition on the normal derivative of the variable varname and on the mesh region re...
 
size_type add_normal_derivative_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname, size_type region)
Adds a normal derivative source term brick :math:F = \int b.
 
size_type add_Kirchhoff_Love_Neumann_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname1, const std::string &dataname2, size_type region)
Adds a Neumann term brick for Kirchhoff-Love model on the variable varname and the mesh region region...
 
size_type add_normal_derivative_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalisation_coeff, size_type region, const std::string &dataname=std::string(), bool R_must_be_derivated=false)
Adds a Dirichlet condition on the normal derivative of the variable varname and on the mesh region re...