Re: [Getfem-users] different mesh_fem definitions for different regions

[Umut Tabak]

Renard Yves wrote:
> Umut Tabak <address@hidden> a écrit :
>
> > Renard Yves wrote:
> > > Dear Umut,
> > >
> > > If your region is a part of a mesh (not a boundary)
Hi again,

I tried to apply a sample analysis where I can assemble the fluid 
part(scalar) and the structural part(displacement). Now I have a problem 
concerning the coupling matrix assembly. I attached the code and the 
necessary msh files. Now the questions are(better to answer these after 
taking a look at the source file attached. In brief, I have 3 mesh 
structures and 3 mesh_fems attached to these mesh structures and the 
operations should be clear for you.)

+ To have this kind of analysis(coupled structural-acoustic), we have to 
set the mesh_fem dimensions to 3 and 1 for the structural and fluid 
parts respectively, if there is only one mesh_fem definition, how can I 
assign different dimensions to different regions of my mesh?

+ The second question follows accordingly, how to assign different 
integration methods to different regions of the mesh? Assembly can be 
done on regions as you suggested before.

+ Without the definition of separate mesh_fem objects, how can I 
allocate memory for the matrices that I am going to use for the assembly 
process? I mean to get the specific dof numbers associated to regions.

+ Lastly, with the attached file, I can assemble structural domain and 
apply boundary conditions on that, also the fluid part, but i have a 
problem related to a term like

\begin{displaymath}
 \bold{K}_{up}=\int_{S}{\bold{N}_{s}}^{T} \bold{n}{\bold{N}_{p}}dS
\end{displaymath}

that I would like to assemble on the interface of the two above defined 
domains.

Shape functions N_s will come from the structural domain and N_p will be 
the contributions of the fluid domain. But i could not assemble this 
term over the surface defined in the mesh, I understand the error 
message but I could not think a way to fix it. The error is ''error: the 
mesh_fem/mesh_im live on different meshes!''.

I would appreciate further help on this.

Thanks and best regards,

Umut

// standard headers
#include <iostream>
#include <vector>
#include <string>
#include <cstdlib>
#include <cassert>
#include <stdexcept>
// gmm++ headers, pay attention that this is a template
// library, so that no build is required.
#include <gmm/gmm.h>
#include <gmm/gmm_inoutput.h>
// getFem++ headers 
#include <getfem/getfem_import.h>
#include <getfem/dal_bit_vector.h>
#include <getfem/getfem_mesh.h>
#include <getfem/getfem_mesh_fem.h>
#include <getfem/getfem_partial_mesh_fem.h>
#include <getfem/getfem_mesh_im.h>
#include <getfem/getfem_integration.h>
#include <getfem/getfem_assembling.h>
#include <getfem/getfem_modeling.h>
#include <getfem/getfem_regular_meshes.h>
//
typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
//using namespace getfem;
//using namespace gmm;
using namespace std;
using bgeot::scalar_type;
//
int main(int argc, char** argv)
{ 
  try{
    // if(argc!=3){
    //   throw std::runtime_error("Error in input parameters\nUsage: <binary> 
<filename.msh> <gmshFormatSpecifier>");
    // }
    getfem::mesh mshGmsh1;
    getfem::mesh mshGmsh2;
    getfem::mesh mshGmsh3;
    //import_mesh(argv[1], argv[2], mshGmsh);
    getfem::import_mesh("gmshv2:./boxSolidDom1.msh", mshGmsh1);
    getfem::import_mesh("gmshv2:./boxSolidDom2.msh", mshGmsh2);
    getfem::import_mesh("gmshv2:./boxSolidDom3.msh", mshGmsh3);
                // optional print out for the convexes(elements) of a specific
                // region
                cout << "Number of elements in region 1 is " << 
mshGmsh1.nb_convex() << endl;
                cout << "Number of elements in region 2 is " << 
mshGmsh2.nb_convex() << endl;
                cout << "Number of elements in region 3 is " << 
mshGmsh3.nb_convex() << endl;
    // can also be an array of vectors
                // dal::bit_vector elemsInRegion1;
                // dal::bit_vector elemsInRegion2;
                // dal::bit_vector elemsInRegion3;
                // element definitions 
                // for (getfem::mr_visitor i(region1); !i.finished(); ++i) 
elemsInRegion1.add(i.cv());
                // for (getfem::mr_visitor i(region2); !i.finished(); ++i) 
elemsInRegion2.add(i.cv());
                // for (getfem::mr_visitor i(region3); !i.finished(); ++i) 
elemsInRegion3.add(i.cv());
    // assign fem definitions to the regions
                getfem::mesh_fem meshFemDom1(mshGmsh1);
                getfem::mesh_fem meshFemDom2(mshGmsh2);
                getfem::mesh_fem meshFemDom3(mshGmsh3);
                // define also a mesh fem for the right hand side 
                // for boundary condition application
                getfem::mesh_fem rhsAssembly(mshGmsh1);
                getfem::mesh_fem lameConsts(mshGmsh1);
                // structural domain 
                meshFemDom1.set_qdim(3);
                
meshFemDom1.set_finite_element(getfem::fem_descriptor("FEM_QK(3, 1)"));
    cout << "Number of dofs of the structural domain: " << meshFemDom1.nb_dof() 
<< endl;
    rhsAssembly.set_qdim(1);
    rhsAssembly.set_finite_element(mshGmsh1.convex_index(), 
                                   getfem::fem_descriptor("FEM_QK(3, 1)"));
    lameConsts.set_qdim(1);
    lameConsts.set_finite_element(mshGmsh1.convex_index(), 
                           getfem::fem_descriptor("FEM_QK(3, 0)"));
                // set integration methods of different domains
    getfem::mesh_im intMethods1(mshGmsh1);
    getfem::pintegration_method ppi1 = 
getfem::int_method_descriptor("IM_HEXAHEDRON(5)");
    intMethods1.set_integration_method(ppi1);
                // define the matrix to use in the assembly process of
    sparse_matrix SM(meshFemDom1.nb_dof(), meshFemDom1.nb_dof());
    sparse_matrix MM(meshFemDom1.nb_dof(), meshFemDom1.nb_dof());
    // define the source term
    vector<scalar_type> B(meshFemDom1.nb_dof());
    // define elastic constants
    vector<double> lambda(lameConsts.nb_dof(), 70e9*0.33/((1+0.33)*(1-2*0.33)));
    vector<double> mu(lameConsts.nb_dof(), 70e9/(2*(1+0.33)));
    // bc values array
    vector<scalar_type> V(meshFemDom1.nb_dof(),0.0);
    V[7] = 10;
    cout << "number of dofs of the RHS term" << rhsAssembly.nb_dof() << endl;
    // assembly process
    std::cout<<" ---  stiffness matrix assembling --- " << std::endl;
    //getfem::asm_stiffness_matrix_for_linear_elasticity(SM, intMethods, 
meshFemDom1, rhsAssembly, lambda, mu);
    getfem::asm_stiffness_matrix_for_linear_elasticity(SM, intMethods1, 
meshFemDom1, lameConsts, lambda, mu);
    std::cout<<" ---  stiffness matrix assembled  --- " << std::endl;
    //
    std::cout<<" ---       mass matrix assembling --- " << std::endl;
    getfem::asm_mass_matrix(MM, intMethods1, meshFemDom1,meshFemDom1);
    std::cout<<" ---       mass matrix assembled  --- " << std::endl;
    //
    std::cout<<" --- source term assembling --- " << std::endl;
    getfem::asm_source_term(B, intMethods1, meshFemDom1, rhsAssembly, V, 2);
    //getfem::asm_source_term(B, intMethods, meshFemDom1, lameConsts, V, 2);
    std::cout<<" ---  source term assembled --- " << std::endl;
    // //
    std::cout<<" --- applying the boundary conditions --- " << std::endl;
    getfem::assembling_Dirichlet_condition(SM, B, meshFemDom1, 2, V);
    getfem::assembling_Dirichlet_condition(MM, B, meshFemDom1, 2, V);
    std::cout<<" --- applied the boundary conditions  --- " << std::endl;
    // dump matrices onto files
                gmm::csc_matrix<double> SM2,MM2;
    string kFileName("KmatDom1.mm");
    string mFileName("MmatDom1.mm");
    gmm::copy(SM, SM2);
    gmm::copy(MM, MM2);
    MatrixMarket_save(mFileName.c_str(), MM2);
    MatrixMarket_save(kFileName.c_str(), SM2);
    // fluid domain, Domain 2, scalar fem, only pressures are considered
                meshFemDom2.set_qdim(1);
    meshFemDom2.set_finite_element(mshGmsh2.convex_index(), 
                                   getfem::fem_descriptor("FEM_QK(3, 1)"));
                // set integration methods of different domains
    getfem::mesh_im intMethods2(mshGmsh2);
    getfem::pintegration_method ppi2 = 
getfem::int_method_descriptor("IM_HEXAHEDRON(5)");
    intMethods2.set_integration_method(ppi2);
    // 
                sparse_matrix SMF(meshFemDom2.nb_dof(), meshFemDom2.nb_dof());
    sparse_matrix MMF(meshFemDom2.nb_dof(), meshFemDom2.nb_dof());
          // assemble fluid only part
    getfem::generic_assembly assem;
    assem.push_mi(intMethods2);
    assem.push_mf(meshFemDom2);
    assem.push_mat(SMF);
    assem.set("M$1(#1,#1)+=sym(comp(Grad(#1).Grad(#1))(:,k,:,k))");
    std::cout<<" ---  Fluid stiffness matrix assembling --- " << std::endl;
    assem.assembly();
    std::cout<<" ---  Fluid stiffness matrix assembled  --- " << std::endl;
    // std::cout<<" ---  stiffness matrix:  " << endl << SM << std::endl;
    std::cout<<" ---  Fluid mass matrix assembling --- " << std::endl;
    getfem::asm_mass_matrix(MMF, intMethods2, meshFemDom2, meshFemDom2);
    std::cout<<" ---  Fluid mass matrix assembled  --- " << std::endl;
    // dump matrices of the fluid part
                gmm::csc_matrix<double> SM3,MM3;
    gmm::copy(SMF, SM3);
    gmm::copy(MMF, MM3);
    string kFileNameFluid("KmatFlui.mm");
    string mFileNameFluid("MmatFlui.mm");
    MatrixMarket_save(mFileNameFluid.c_str(), MM3);
    MatrixMarket_save(kFileNameFluid.c_str(), SM3);
                // set integration methods for the coupling surface 
    getfem::mesh_im intMethods3(mshGmsh3);
    getfem::pintegration_method ppi3 = 
getfem::int_method_descriptor("IM_QUAD(3)");
    intMethods3.set_integration_method(ppi3);
    // allocate the coupling matrix
                sparse_matrix SMCoupling(gmm::mat_nrows(SM), 
gmm::mat_nrows(SMF));
    //  assemble coupling terms
    getfem::generic_assembly assemCoupling;
    assemCoupling.push_mi(intMethods3);
                assemCoupling.push_mf(meshFemDom1); // structural domain #1
    assemCoupling.push_mf(meshFemDom2); // fluid domain      #2
    assemCoupling.push_mat(SMCoupling); 
    
assemCoupling.set("M$1(#1,#2)+=sym(comp(Base(#1).Normal().Base(#2))(:,k,:,k))");
    std::cout<<" ---  Coupling matrix assembling --- " << std::endl;
    assemCoupling.assembly();
    std::cout<<" ---  Coupling matrix assembled  --- " << std::endl;
    // dump matrices of the fluid part
                gmm::csc_matrix<double> SMC;
    gmm::copy(SMCoupling, SMC);
    string kFileNameCoupling("KmatC.mm");
    MatrixMarket_save(kFileNameCoupling.c_str(),SMC);
  }
  catch(std::exception &e){
    std::cout << e.what() << std::endl;
  }
  return EXIT_SUCCESS;
}

boxSolidDom1.msh
Description: Mesh model

boxSolidDom2.msh
Description: Mesh model

boxSolidDom3.msh
Description: Mesh model