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

[Umut Tabak]

Renard Yves wrote:
> Dear Umut,
>
> I your exemple, you should define two mesh_fems as follows :
>
>
>
>    getfem::mesh_fem meshFemCore1(mshGmsh);
>    // region 1, scalar domain, pressure field, dimension 1
>    for(getfem::mr_visitor i(region1); !i.finished(); ++i){
>        
> meshFemCore1.set_finite_element(i.cv(),getfem::fem_descriptor("FEM_QK(3, 
> 1)"));
>    }
>    meshFemCore.set_qdim(1); // not really necessary.
>
>    getfem::mesh_fem meshFemCore2(mshGmsh);
>    // region 2, vectorial domain, displacement field, dimension 3
>    for(getfem::mr_visitor i(region2); !i.finished(); ++i){
>        
> meshFemCore2.set_finite_element(i.cv(),getfem::fem_descriptor("FEM_QK(3, 
> 1)"));
>    }
>    meshFemCore.set_qdim(3);
>
>
>
> Then, you can use these two mesh_fem objects in the assembly procedure 
> in order to compute your coupling term.
Dear Professor Renard,

I guess with your great help, I am close to what I would like to get, I 
updated my code and can define regions on one mesh. The last problem I 
am having is that I can not assemble, integral term on the interface 
surface, well not really, I can run the code without problems but the 
point is that I do not get any terms. For this integral, i need to 
integrate that term over quad regions, but using the predefined 
mesh_fems, this does not seem possible to me, or I am still having some 
mistakes, the code is attached, region 3 is the interface region inside 
the mesh domain on which I would like to assemble my coupling terms. Am 
I experiencing a problem due to the normal definition, or is there a 
problem in the index expressions?

Best regards,
Umut

boxSolid.msh
Description: Mesh model

// 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 mshGmsh;
    //import_mesh(argv[1], argv[2], mshGmsh);
    // there is one mesh definition and different mesh_fem definitions
    // for different regions in definition
    getfem::import_mesh("gmshv2:./boxSolid.msh", mshGmsh);
    if(mshGmsh.has_region(1) && mshGmsh.has_region(2) && mshGmsh.has_region(3) 
&& mshGmsh.has_region(4) ){
                        std::cout << "Mesh has the specified 4 region 
definitions" << std::endl;
                }
                // 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(mshGmsh);
    // set definitions for region 1
    for(getfem::mr_visitor i(mshGmsh.region(1)); !i.finished(); ++i){
                        
meshFemDom1.set_finite_element(i.cv(),getfem::fem_descriptor("FEM_QK(3, 1)"));
    }
    // scalar field dimension is set to 1.
    meshFemDom1.set_qdim(1);
    // set defintions for region 2
    getfem::mesh_fem meshFemDom2(mshGmsh);
    for(getfem::mr_visitor i(mshGmsh.region(2)); !i.finished(); ++i){
                        
meshFemDom2.set_finite_element(i.cv(),getfem::fem_descriptor("FEM_QK(3, 1)"));
    }
    // vectorial domain so that the dimension is set to 3.
    meshFemDom2.set_qdim(3);
    // for boundary condition application
    getfem::mesh_fem rhsAssembly(mshGmsh);
    getfem::mesh_fem lameConsts(mshGmsh);
    // structural domain 
    cout << "Number of dofs of the structural domain: " << meshFemDom2.nb_dof() 
<< endl;
    //
                for(getfem::mr_visitor i(mshGmsh.region(2)); !i.finished(); 
++i){
                        
rhsAssembly.set_finite_element(i.cv(),getfem::fem_descriptor("FEM_QK(3, 1)"));
                }
    rhsAssembly.set_qdim(1);
    //
                for(getfem::mr_visitor i(mshGmsh.region(2)); !i.finished(); 
++i){
                        
lameConsts.set_finite_element(i.cv(),getfem::fem_descriptor("FEM_QK(3, 0)"));
                }
                lameConsts.set_qdim(1);
                // set integration methods of different domains
    getfem::mesh_im intMethod(mshGmsh);
    getfem::pintegration_method ppi1 = 
getfem::int_method_descriptor("IM_HEXAHEDRON(5)");
    intMethod.set_integration_method(ppi1);
                // define the matrix to use in the assembly process of
    sparse_matrix SM(meshFemDom2.nb_dof(), meshFemDom2.nb_dof());
    sparse_matrix MM(meshFemDom2.nb_dof(), meshFemDom2.nb_dof());
    // define the source term
    vector<scalar_type> B(meshFemDom2.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(meshFemDom2.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, intMethod, 
meshFemDom2, lameConsts, lambda, mu);
    std::cout<<" ---  stiffness matrix assembled  --- " << std::endl;
    //
    std::cout<<" ---       mass matrix assembling --- " << std::endl;
    getfem::asm_mass_matrix(MM, intMethod, meshFemDom2,meshFemDom2);
    std::cout<<" ---       mass matrix assembled  --- " << std::endl;
    //
    std::cout<<" --- source term assembling --- " << std::endl;
    getfem::asm_source_term(B, intMethod, meshFemDom2, rhsAssembly, V);
    //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, meshFemDom2, 4, V);
    getfem::assembling_Dirichlet_condition(MM, B, meshFemDom2, 4, V);
    std::cout<<" --- applied the boundary conditions  --- " << std::endl;
    // dump matrices onto files
                gmm::csc_matrix<double> SM2,MM2;
    string kFileName("KmatDom2.mm");
    string mFileName("MmatDom2.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
                sparse_matrix SMF(meshFemDom1.nb_dof(), meshFemDom1.nb_dof());
    sparse_matrix MMF(meshFemDom1.nb_dof(), meshFemDom1.nb_dof());
          // assemble fluid only part
    getfem::generic_assembly assem;
    assem.push_mi(intMethod);
    assem.push_mf(meshFemDom1);
    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, intMethod, meshFemDom1, meshFemDom1);
    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);
                std::cout << "Matrices save in matrix market format" << 
std::endl;
    // 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(intMethod);
                assemCoupling.push_mf(meshFemDom2); // structural domain #1
    assemCoupling.push_mf(meshFemDom1); // fluid domain      #2
    assemCoupling.push_mat(SMCoupling); 
    assemCoupling.set("M$1(#1,#2)+=comp(vBase(#1).Normal().Base(#2))(:,i,i,:)");
    std::cout << " ---  Coupling matrix assembling --- " << std::endl;
    assemCoupling.assembly(mshGmsh.region(3));
    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);
    std::cout<<" ---  Coupling matrix saved  --- " << std::endl;
  }
  catch(std::exception &e){
    std::cout << e.what() << std::endl;
  }
  return EXIT_SUCCESS;
}