[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Getfem-users] different mesh_fem definitions for different regions
From: |
Umut Tabak |
Subject: |
Re: [Getfem-users] different mesh_fem definitions for different regions |
Date: |
Thu, 30 Jul 2009 14:09:18 +0200 |
User-agent: |
Thunderbird 2.0.0.21 (X11/20090409) |
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