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...