Welcome to the NetCologne GmbH open source mirroring service!
This machine mirrors various open-source projects.
20 Gbit/s uplink.
If there are any issues or you want another project mirrored, please contact
mirror-service -=AT=- netcologne DOT de !
GetFEM: src/getfem/getfem_level_set_contact.h Source File
Go to the documentation of this file.
47 #ifndef GETFEM_LEVEL_SET_CONTACT_H__
48 #define GETFEM_LEVEL_SET_CONTACT_H__
56 namespace level_set_contact {
63 using getfem::scalar_type;
64 using getfem::model_real_plain_vector;
65 typedef getfem::model_real_plain_vector plain_vector;
66 typedef getfem::model_real_sparse_matrix sparse_matrix;
82 getfem::model_real_plain_vector RHS(mf.
nb_dof());
85 getfem::base_node Pmin(
mesh.dim()),Pmax(
mesh.dim()),range(
mesh.dim());
87 gmm::add(Pmax,gmm::scaled(Pmin,-1.0),range);
88 getfem::scalar_type mesh_size = *(std::max_element(range.begin(),range.end()));
89 gmm::fill(RHS,mesh_size*5.0);
106 GMM_TRACE2(
"building scalar level set with laplace equation..");
114 GMM_TRACE2(
"..done");
121 const std::string var_name;
133 inline std::string get_var_name()
const {
return var_name;}
134 inline mesh& get_mesh() {
return own_mesh;}
135 inline const mesh& get_mesh()
const {
return own_mesh;}
136 inline const mesh_fem& get_mesh_fem()
const {
return own_mesh_fem;}
137 inline const model& get_model()
const {
return md;}
138 inline bool is_mesh_deformed()
const {
return is_deformed;}
159 std::string _ls_name);
160 inline std::string get_ls_name()
const {
return ls_name;}
161 inline const plain_vector& ls_values()
const
163 inline plain_vector& ls_values()
166 template<
class VECTOR>
void set_level_set(
const VECTOR& ls)
183 const std::string mult_name;
187 dal::bit_vector old_contact_elm_list;
188 dal::bit_vector pre_old_ct_list;
190 mutable std::shared_ptr<mesh_im> pmim_contact;
192 mutable std::shared_ptr<mesh_fem> pinterpolated_fem;
193 mutable std::shared_ptr<mesh_fem> pinterpolated_fem_U;
194 mutable std::shared_ptr<gmm::unsorted_sub_index> slave_ls_dofs;
195 mutable std::shared_ptr<gmm::unsorted_sub_index> slave_U_dofs;
199 mutable bool members_are_computed;
200 mutable bool init_cont_detect_done;
204 inline const mesh_fem& slave_scalar_fem()
const {
205 if (dependecies_changed())
update();
206 return *pinterpolated_fem;
208 inline const mesh_fem& slave_vector_fem()
const {
209 if (dependecies_changed())
update();
210 return *pinterpolated_fem_U;
212 inline const gmm::unsorted_sub_index& slave_scalar_dofs()
const {
213 if (dependecies_changed())
update();
214 return *slave_ls_dofs;
216 inline const gmm::unsorted_sub_index& slave_vector_dofs()
const {
217 if (dependecies_changed())
update();
218 return *slave_U_dofs;
220 inline const mesh_im& contact_mesh_im()
const {
221 if (dependecies_changed())
update();
222 return *pmim_contact;
226 {
return ACTIVE_CONTACT_REGION;}
228 inline const std::string& get_mult_name()
const
231 inline size_type num_of_integr_elems()
const {
return n_integrated_elems;}
233 inline bool dependecies_changed()
const
234 {
return !members_are_computed;}
235 inline void force_update()
const
236 {members_are_computed=
false;}
291 const std::string mult_int_method;
292 size_type BOUNDARY_ELEMENTS, VOLUME_ELEMENTS;
293 std::vector<size_type> face_to_belem_ind;
294 static std::vector<master_contact_body*> masters;
295 std::map<std::string, std::shared_ptr<contact_pair_info> > contact_table;
296 std::map<size_type,face_type> border_faces;
308 enum contact_integration{PER_ELEMENT=1,REGULARIZED_LEVEL_SET=2};
340 const std::string& _var_name,
351 const std::string& _var_name,
353 const std::string& _mult_int_method,
354 contact_integration _integration = PER_ELEMENT,
355 scalar_type _regularized_tollerance = 1e-6,
356 scalar_type _small_weight_multiplier = 0.001,
357 scalar_type _max_contact_angle = 45);
368 GMM_ASSERT1(mult_mim_order!=
size_type(-1),
369 "master body was not created with "
370 "order of integration for contact area");
371 return mult_mim_order;
379 GMM_ASSERT1(mult_mim_order==
size_type(-1),
380 "master body was not created with integration "
381 "method for contact area");
387 {
return VOLUME_ELEMENTS;}
392 {
return BOUNDARY_ELEMENTS;}
415 inline void update_for_slave(std::string slave_var_name)
416 { contact_table[slave_var_name]->update(); }
433 enum update_depth{DEFORM_MESHES_ONLY,FULL_UPDATE};
438 std::shared_ptr<getfem::temporary_mesh_deformator<> > def_master;
439 std::shared_ptr<getfem::temporary_mesh_deformator<> > def_slave;
445 update_depth ud = FULL_UPDATE);
490 const model::varnamelist &vl,
491 const model::varnamelist &dl,
492 const model::mimlist &mims,
493 model::real_matlist &matl,
494 model::real_veclist &vecl,
495 model::real_veclist &,
497 build_version version)
const;
514 bgeot::multi_index sizes_;
524 dim(_mcb.get_mesh().dim()) {
526 GMM_ASSERT1(dim==2 || dim==3,
"NormalTerm: wrong space dimension ");
527 GMM_ASSERT1(version==1 || version==2,
"NormalTerm:: wrong version ");
544 const bgeot::multi_index &sizes(
size_type)
const {
return sizes_;};
560 const plain_vector &LS_U;
561 scalar_type m_Epsilon;
563 bgeot::multi_index sizes_;
569 const plain_vector &LS_U_,
570 scalar_type epsilon=1e-9,
571 scalar_type small_h_=0);
572 const bgeot::multi_index &sizes(
size_type)
const;
575 scalar_type hRegularized(scalar_type x, scalar_type epsion, scalar_type small);
583 bgeot::multi_index sizes_;
587 const bgeot::multi_index &sizes(
size_type)
const;
594 template<
typename MAT,
typename VECT>
595 void asm_level_set_contact_tangent_matrix(
596 std::vector<MAT>& matl,
597 const master_contact_body& mcb,
598 const slave_contact_body& scb,
610 const std::string& mult_name =
611 mcb.get_pair_info(scb.get_var_name()).get_mult_name();
612 const std::string ls_name =
"ls_on"+mcb.get_var_name()+
"_from_"+scb.get_var_name();
615 const mesh_fem& mf_U_line = mcb.get_mesh_fem();
616 const mesh_fem& mf_lambda = mcb.get_model().mesh_fem_of_variable(mult_name);
617 const mesh_fem& mf_interpolate =
618 mcb.get_pair_info(scb.get_var_name()).slave_scalar_fem();
619 const mesh_fem& mf_U_interpolate =
620 mcb.get_pair_info(scb.get_var_name()).slave_vector_fem();
621 const mesh_fem& mf_master_ls = mcb.get_model().mesh_fem_of_variable(ls_name);
622 const mesh_im& mim_line =
623 mcb.get_pair_info(scb.get_var_name()).contact_mesh_im();
626 plain_vector LS_small(mf_interpolate.nb_dof());
627 gmm::copy(gmm::sub_vector(scb.ls_values(),
628 mcb.get_pair_info(scb.get_var_name()).slave_scalar_dofs()),LS_small);
631 NormalTerm R_matrix(mcb,2);
634 std::shared_ptr<getfem::nonlinear_elem_term> integration;
635 if (mcb.integration==master_contact_body::REGULARIZED_LEVEL_SET) {
636 integration = std::make_shared<HFunction>
637 (mf_master_ls,mcb.get_model().real_variable(ls_name),
638 mcb.regularized_tollerance,mcb.small_weight_multiplier);
639 }
else { integration = std::make_shared<Unity>(mf_master_ls); }
643 sparse_matrix Kms_small(mf_U_interpolate.nb_dof(), mf_U_line.nb_dof());
644 sparse_matrix Kss_small(mf_U_interpolate.nb_dof(),mf_U_interpolate.nb_dof());
645 sparse_matrix Ksl_small(mf_U_interpolate.nb_dof(),mf_lambda.nb_dof());
653 "Kmm1 = comp(Base(#1).Grad(#3).vBase(#2).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,k,:,k,m,n,:,n,m,1).L(i).F(j);"
654 "Kmm2 = comp(Base(#1).NonLin$1(#2).vGrad(#2).Grad(#3).vBase(#2).NonLin$2(#5))(i,m,n,:,n,m,j,k,:,k,1).L(i).F(j);"
655 "Kmm3 = comp(Base(#1).Base(#3).NonLin$1(#2).vGrad(#2).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,k,l,:,l,k,m,n,:,n,m,1).L(i).F(j);"
656 "Kmm4 = (-1.0)*comp(Base(#1).Base(#3).NonLin$1(#2).vGrad(#2).vGrad(#2).NonLin$2(#5))(i,j,m,n,:,n,l,:,l,m,1).L(i).F(j);"
657 "M$1(#2,#2)+= sym(Kmm1+Kmm2+Kmm3+Kmm4);"
658 "Ksm1=(-1.0)*comp(Base(#1).Grad(#3).vBase(#4).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,k,:,k,m,n,:,n,m,1).L(i).F(j);"
659 "Ksm2=(-1.0)*comp(Base(#1).Grad(#3).vGrad(#4).vBase(#2).NonLin$2(#5))(i,j,m,:,m,n,:,n,1).L(i).F(j);"
660 "M$2(#4,#2)+= Ksm1+Ksm2;"
661 "Kml1=comp(Base(#3).NonLin$1(#2).vGrad(#2).Base(#1).NonLin$2(#5))(i,m,n,:,n,m,:,1).F(i);"
662 "Kml2=comp(Grad(#3).vBase(#2).Base(#1).NonLin$2(#5))(i,j,:,j,:,1).F(i);"
663 "M$3(#2,#1)+= Kml1+Kml2;"
664 "Kss_part = comp(Base(#1).Grad(#3).vGrad(#4).vBase(#4).NonLin$2(#5))(i,j,m,:,m,n,:,n,1).L(i).F(j);"
665 "M$4(#4,#4)+=sym(Kss_part{1,2}+Kss_part{2,1});"
666 "M$5(#4,#1)+=(-1.0)*comp(Grad(#3).vBase(#4).Base(#1).NonLin$2(#5))(i,k,:,k,:,1).F(i);"
670 assem_boundary.
push_mi(mim_line);
671 assem_boundary.
push_mf(mf_lambda);
672 assem_boundary.
push_mf(mf_U_line);
673 assem_boundary.
push_mf(mf_interpolate);
674 assem_boundary.
push_mf(mf_U_interpolate);
675 assem_boundary.
push_mf(mf_master_ls);
685 assem_boundary.
assembly(contact_region);
688 const gmm::sub_interval& Um_dof = gmm::sub_interval(0,mf_U_line.nb_dof());
689 const gmm::unsorted_sub_index& Us_dof =
690 mcb.get_pair_info(scb.get_var_name()).slave_vector_dofs();
691 const gmm::sub_interval& LM_dof = gmm::sub_interval(0,mf_lambda.nb_dof());
692 gmm::copy(gmm::transposed(Kms_small),gmm::sub_matrix(Kms,Um_dof,Us_dof));
693 gmm::copy(Kss_small,gmm::sub_matrix(Kss,Us_dof,Us_dof));
694 gmm::copy(Ksl_small,gmm::sub_matrix(Ksl,Us_dof,LM_dof));
698 template<
typename VECT0,
typename VECT1>
699 void asm_level_set_contact_rhs(
700 std::vector<VECT0>& vecl,
701 const master_contact_body& mcb,
702 const slave_contact_body& scb,
707 VECT0& RHS_Um = vecl[0];
708 VECT0& RHS_Us = vecl[1];
709 VECT0& RHS_LM = vecl[2];
713 const std::string& mult_name =
714 mcb.get_pair_info(scb.get_var_name()).get_mult_name();
715 const std::string ls_name =
"ls_on"+mcb.get_var_name()+
"_from_"+scb.get_var_name();
718 const mesh_fem& mf_U_line = mcb.get_mesh_fem();
719 const mesh_fem& mf_lambda =
720 mcb.get_model().mesh_fem_of_variable(mult_name);
721 const mesh_fem& mf_interpolate =
722 mcb.get_pair_info(scb.get_var_name()).slave_scalar_fem();
723 const mesh_fem& mf_U_interpolate =
724 mcb.get_pair_info(scb.get_var_name()).slave_vector_fem();
725 const mesh_fem& mf_master_ls = mcb.get_model().mesh_fem_of_variable(ls_name);
726 const mesh_im& mim_line =
727 mcb.get_pair_info(scb.get_var_name()).contact_mesh_im();
730 plain_vector LS_small(mf_interpolate.nb_dof());
731 gmm::copy(gmm::sub_vector(scb.ls_values(),
732 mcb.get_pair_info(scb.get_var_name()).slave_scalar_dofs()),LS_small);
735 NormalTerm R_matrix(mcb,2);
738 std::shared_ptr<getfem::nonlinear_elem_term> integration;
739 if (mcb.integration==master_contact_body::REGULARIZED_LEVEL_SET) {
740 integration = std::make_shared<HFunction>
741 (mf_master_ls,mcb.get_model().real_variable(ls_name),
742 mcb.regularized_tollerance,mcb.small_weight_multiplier);
743 }
else { integration = std::make_shared<Unity>(mf_master_ls); }
746 plain_vector RHS_Us_small(mf_U_interpolate.nb_dof());
752 "RHS_L_Us_1=comp(Base(#1).Base(#3).NonLin$1(#2).vGrad(#2).NonLin$2(#5))(i,j,m,n,:,n,m,1).L(i).F(j);"
753 "RHS_L_Us_2=comp(Base(#1).Grad(#3).vBase(#2).NonLin$2(#5))(i,j,k,:,k,1).L(i).F(j);"
754 "V$1(#2)+=RHS_L_Us_1+RHS_L_Us_2;"
755 "V$2(#4)+=(-1.0)*comp(Base(#1).Grad(#3).vBase(#4).NonLin$2(#5))(i,j,k,:,k,1).L(i).F(j);"
756 "V$3(#1)+=comp(Base(#1).Base(#3).NonLin$2(#5))(:,i,1).F(i);"
758 assem_boundary.
push_mi(mim_line);
759 assem_boundary.
push_mf(mf_lambda);
760 assem_boundary.
push_mf(mf_U_line);
761 assem_boundary.
push_mf(mf_interpolate);
762 assem_boundary.
push_mf(mf_U_interpolate);
763 assem_boundary.
push_mf(mf_master_ls);
769 assem_boundary.
push_vec(RHS_Us_small);
771 assem_boundary.
assembly(contact_region);
774 const gmm::unsorted_sub_index& Us_dof =
775 mcb.get_pair_info(scb.get_var_name()).slave_vector_dofs();
776 gmm::copy(RHS_Us_small, gmm::sub_vector(RHS_Us,Us_dof));
782 typedef void(*SOLVE_FUNCTION)(
785 getfem::rmodel_plsolver_type solver,
786 getfem::abstract_newton_line_search &ls);
804 const std::string& lsolver,
805 getfem::abstract_newton_line_search &ls);
809 #endif //GETFEM_LEVEL_SET_CONTACT_H__
Contact body that will be projected on the boundary of the master.
Standard solvers for model bricks.
const contact_pair_info & get_pair_info(const std::string &slave_var_name) const
access to a structure that contains all the info about contact pair between this master and a slave,...
master_contact_body(model &_md, const std::string &_var_name, size_type _mult_order, size_type _mult_mim_order)
create master contact body with a model, name where masters displacements are defined,...
const mesh_region region(size_type id) const
Return the region of index 'id'.
bool master_contact_changed(void)
contact detection for all slaves
getfem::pintegration_method contact_int_method() const
integration method on the contact surface, use it when the master is created with a specific integrat...
void clear_contact_history(void)
clearing previous contact elem lists
size_t size_type
used as the common size type in the library
pintegration_method int_method_descriptor(std::string name, bool throw_if_not_found=true)
Get an integration method from its name .
Describe an integration method linked to a mesh.
slave_contact_body(model &_md, const std::string &_var_name, mesh_im *_pmim)
default constructor.
void add_initialized_fem_data(const std::string &name, const mesh_fem &mf, const VECT &v)
Add an initialized fixed size data to the model, assumed to be a vector field if the size of the vect...
void add_slave(slave_contact_body &scb, size_type slave_contact_region=-1)
associate a slave contact body with this master.
void push_mf(const mesh_fem &mf_)
Add a new mesh_fem.
void push_mat(const MAT &m)
Add a new output matrix (fake const version..)
void sup_region(size_type b)
Remove the region of index b.
Definition of basic exceptions.
Describe a finite element method linked to a mesh.
static size_type free_region_id(const getfem::mesh &m)
Extract the next region number that does not yet exists in the mesh.
`‘Model’' variables store the variables, the data and the description of a model.
gmm::uint16_type short_type
used as the common short type integer in the library
void push_data(const VEC &d)
Add a new data (dense array)
The Iteration object calculates whether the solution has reached the desired accuracy,...
structure passed as the argument of fem interpolation functions.
void assembly(const mesh_region ®ion=mesh_region::all_convexes())
do the assembly on the specified region (boundary or set of convexes)
static bool any_contact_change()
contact detection for all masters/slave couples
Model representation in Getfem.
"iterator" class for regions.
void add_fem_variable(const std::string &name, const mesh_fem &mf, size_type niter=1)
Add a variable being the dofs of a finite element method to the model.
const contact_integration integration
integration approach for contact elements that are partially crossed by level sets: PER_ELEMENT - a w...
Generic assembly of vectors, matrices.
const scalar_type small_weight_multiplier
in case of REGULARIZED_LEVEL_SET this value scales weight of Gauss points that have negative level se...
void bounding_box(base_node &Pmin, base_node &Pmax) const
Return the bounding box [Pmin - Pmax] of the mesh.
const scalar_type regularized_tollerance
width of transition for a regularazied Heaviside function in case of REGULARIZED_LEVEL_SET
size_type contact_mim_order() const
order of integration of boundary contact terms
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
const mesh_fem & mesh_fem_of_variable(const std::string &name) const
Gives the access to the mesh_fem of a variable if any.
std::shared_ptr< mesh_im > build_mesh_im_on_boundary(size_type region)
return a pointer to mesh_im used for contact surface calculations
base class for the master and the slave contact bodies.
void push_vec(VEC &v)
Add a new output vector.
void push_mi(const mesh_im &im_)
Add a new mesh_im.
size_type APIDECL add_Dirichlet_condition_with_penalization(model &md, const mesh_im &mim, const std::string &varname, scalar_type penalization_coeff, size_type region, const std::string &dataname=std::string(), const mesh_fem *mf_mult=0)
Add a Dirichlet condition on the variable varname and the mesh region region.
void offset_level_set(scalar_type off)
adds a fixed value "off" to the level set field
Describe a mesh (collection of convexes (elements) and points).
abstract class for integration of non-linear terms into the mat_elem computations the nonlinear term ...
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
The virtual brick has to be derived to describe real model bricks.
structure used to hold a set of convexes and/or convex faces.
face_type ext_face_of_elem(size_type cv) const
gives a face, corresponding to newly created boundary element
size_type APIDECL add_source_term_brick(model &md, const mesh_im &mim, const std::string &varname, const std::string &dataexpr, size_type region=size_type(-1), const std::string &directdataname=std::string())
Add a source term on the variable varname.
size_type volume_region() const
region of all volume elements without the boundary
void APIDECL outer_faces_of_mesh(const mesh &m, const dal::bit_vector &cvlst, convex_face_ct &flist)
returns a list of "exterior" faces of a mesh (i.e.
Master contact body which surface will be used to project contact stresses and stiffness terms.
void push_nonlinear_term(pnonlinear_elem_term net)
Add a new non-linear term.
void standard_solve(model &md, gmm::iteration &iter, rmodel_plsolver_type lsolver, abstract_newton_line_search &ls)
A default solver for the model brick system.
size_type APIDECL add_Laplacian_brick(model &md, const mesh_im &mim, const std::string &varname, size_type region=size_type(-1))
Add a Laplacian term on the variable varname (in fact with a minus : :math:-\text{div}(\nabla u)).
const scalar_type max_contact_angle
if the angle (in degrees) between contact element and level set contour exceed this value,...
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
model_real_plain_vector & set_real_variable(const std::string &name, size_type niter) const
Gives the write access to the vector value of a variable.
size_type boundary_region() const
boundary elements, added after creation of the master contact body
const model_real_plain_vector & real_variable(const std::string &name, size_type niter) const
Gives the access to the vector value of a variable.
Miscelleanous assembly routines for common terms. Use the low-level generic assembly....
static void clear_all_contact_history()
should be used in the beginning of a step to clean data structures that store previous contact elemen...
const size_type mult_mf_order
approximation order for Lagrange multiplier on the contact surface