29 using namespace MoFEM;
32 #include <petsctime.h>
34 #include <boost/numeric/ublas/vector_proxy.hpp>
46 #error "MoFEM need to be compiled with ADOL-C"
49 #include <adolc/adolc.h>
52 #include <boost/numeric/ublas/symmetric.hpp>
56 #include <boost/ptr_container/ptr_map.hpp>
62 double circumcenter[3],
double *xi,
double *
eta,
double *zeta);
64 double circumcenter[3],
double *xi,
double *
eta);
69 #include <moab/AdaptiveKDTree.hpp>
77 static char help[] =
"...\n\n";
100 double xcoord, ycoord, zcoord;
103 ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
112 typedef multi_index_container<
116 tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
119 tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
122 tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
125 tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
128 tag<Composite_xyzcoord>,
131 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
132 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
133 member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
136 int main(
int argc,
char *argv[]) {
145 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
148 PetscBool flg = PETSC_TRUE;
151 if(flg != PETSC_TRUE) {
152 SETERRQ(PETSC_COMM_SELF,
MOFEM_NOT_FOUND,
"*** ERROR -my_file (MESH FILE NEEDED)");
157 if(flg != PETSC_TRUE) {
165 if(flg != PETSC_TRUE) {
173 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
174 if(pcomm == NULL) pcomm =
new ParallelComm(&moab,PETSC_COMM_WORLD);
177 PetscLogDouble t1,t2;
178 PetscLogDouble v1,v2;
179 ierr = PetscTime(&v1); CHKERRQ(
ierr);
180 ierr = PetscGetCPUTime(&t1); CHKERRQ(
ierr);
186 vector<BitRefLevel> bit_levels;
189 int def_meshset_info[2] = {0,0};
190 rval = moab.tag_get_handle(
191 "MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info
196 if(meshset_data[0]==0) {
201 bit_levels.push_back(
BitRefLevel().set(meshset_data[0]-1));
210 Range preriodic_prisms;
218 ierr = m_field.get_cubit_msId_entities_by_dimension(101,
SIDESET,2,SurTrisNeg,
true); CHKERRQ(
ierr);
219 ierr = PetscPrintf(PETSC_COMM_WORLD,
"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(
ierr);
221 double TriCen[3], coords_Tri[9];
224 double roundfact=1e3;
226 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
235 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
236 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
237 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
240 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
241 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
242 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
246 Face_CenPos_Handle_varNeg.insert(
Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
253 ierr = m_field.get_cubit_msId_entities_by_dimension(102,
SIDESET,2,SurTrisPos,
true); CHKERRQ(
ierr);
254 ierr = PetscPrintf(PETSC_COMM_WORLD,
"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(
ierr);
255 for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
264 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
265 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
266 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
269 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
270 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
271 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
275 Face_CenPos_Handle_varPos.insert(
Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
279 double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
280 typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
281 typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
282 typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
283 Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
284 Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
285 Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
288 XcoordMin_it=Face_CenPos_Handle_varNeg.get<
xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
289 XcoordMax_it=Face_CenPos_Handle_varPos.get<
xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
290 YcoordMin_it=Face_CenPos_Handle_varNeg.get<
ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
291 YcoordMax_it=Face_CenPos_Handle_varPos.get<
ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
292 ZcoordMin_it=Face_CenPos_Handle_varNeg.get<
zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
293 ZcoordMax_it=Face_CenPos_Handle_varPos.get<
zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
301 typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
302 Tri_Hand_iterator Tri_Neg;
303 typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
304 xyzcoord_iterator Tri_Pos;
305 double XPos, YPos, ZPos;
309 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
311 Tri_Neg=Face_CenPos_Handle_varNeg.get<
Tri_Hand_tag>().find(*it);
315 if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
316 if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
317 if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
318 Tri_Pos=Face_CenPos_Handle_varPos.get<
Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
325 vector<EntityHandle> TriNodesNeg, TriNodesPos;
326 double CoordNodeNeg[9], CoordNodePos[9];
331 for(
int ii=0; ii<3; ii++){
332 PrismNodes[ii]=TriNodesNeg[ii];
342 double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
343 for(
int ii=0; ii<3; ii++){
344 if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
345 if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
346 if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
347 for(
int jj=0; jj<3; jj++){
349 XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
350 YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
351 ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
354 XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
355 XNodePos=round(XNodePos*roundfact)/roundfact;
356 YNodePos=round(YNodePos*roundfact)/roundfact;
357 ZNodePos=round(ZNodePos*roundfact)/roundfact;
359 if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
360 PrismNodes[3+ii]=TriNodesPos[jj];
366 double CoordNodesPrisms[18];
375 preriodic_prisms.insert(PeriodicPrism);
387 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),
BitRefLevel().set(),out_meshset); CHKERRQ(
ierr);
388 Range LatestRefinedTets;
389 rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,
true);
CHKERRQ_MOAB(
rval);
390 Range LatestRefinedPrisms;
391 rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,
true);
CHKERRQ_MOAB(
rval);
393 Range prims_on_problem_bit_level;
394 ierr = m_field.get_entities_by_type_and_ref_level(
395 problem_bit_level,
BitRefLevel().set(),MBPRISM,prims_on_problem_bit_level
400 rval = moab.add_entities(meshset_prims_on_problem_bit_level,prims_on_problem_bit_level);
CHKERRQ_MOAB(
rval);
401 ierr = m_field.seed_ref_level_MESHSET(meshset_prims_on_problem_bit_level,
BitRefLevel().set()); CHKERRQ(
ierr);
414 if(it->getName().compare(0,12,
"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
429 const double coords[] = {0,0,0};
431 Range range_no_field_vertex;
432 range_no_field_vertex.insert(no_field_vertex);
438 Range surface_negative;
439 ierr = m_field.get_cubit_msId_entities_by_dimension(101,
SIDESET,2,surface_negative,
true); CHKERRQ(
ierr);
440 ierr = PetscPrintf(PETSC_COMM_WORLD,
"number of SideSet 101 = %d\n",surface_negative.size()); CHKERRQ(
ierr);
452 const double coords[] = {0,0,0};
454 Range range_no_field_vertex;
455 range_no_field_vertex.insert(no_field_vertex);
479 PetscBool fo_boundary = PETSC_FALSE;
483 if(it->getName().compare(0,12,
"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
487 rval = moab.get_adjacencies(tris,1,
false,tris_edges,moab::Interface::UNION);
CHKERRQ_MOAB(
rval);
510 ierr = m_field.
loop_dofs(
"MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(
ierr);
512 boost::shared_ptr<Hooke<adouble> > hooke_adouble(
new Hooke<adouble>());
513 boost::shared_ptr<Hooke<double> > hooke_double(
new Hooke<double>());
517 if(it->getName() !=
"MAT_ELASTIC_1")
continue;
519 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(
ierr);
520 int id = it->getMeshsetId();
523 meshset,MBTET,iso_elastic.
setOfBlocks[
id].tEts,
true
528 iso_elastic.
setOfBlocks[id].materialDoublePtr = hooke_double;
529 iso_elastic.
setOfBlocks[id].materialAdoublePtr = hooke_adouble;
533 ierr = iso_elastic.
setOperators(
"DISPLACEMENT",
"MESH_NODE_POSITIONS",
false,
true); CHKERRQ(
ierr);
541 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
542 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
543 bool trans_iso_blocks =
false;
546 string name = it->getName();
547 if (name.compare(0,20,
"MAT_ELASTIC_TRANSISO") == 0) {
548 trans_iso_blocks =
true;
549 int id = it->getMeshsetId();
551 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(
ierr);
552 tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
553 tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
555 tranversly_isotropic_adouble_ptr_map.at(
id)->E_p = mydata.
data.Youngp;
556 tranversly_isotropic_double_ptr_map.at(
id)->E_p = mydata.
data.Youngp;
557 tranversly_isotropic_adouble_ptr_map.at(
id)->E_z = mydata.
data.Youngz;
558 tranversly_isotropic_double_ptr_map.at(
id)->E_z = mydata.
data.Youngz;
559 tranversly_isotropic_adouble_ptr_map.at(
id)->nu_p = mydata.
data.Poissonp;
560 tranversly_isotropic_double_ptr_map.at(
id)->nu_p = mydata.
data.Poissonp;
561 tranversly_isotropic_adouble_ptr_map.at(
id)->nu_pz = mydata.
data.Poissonpz;
562 tranversly_isotropic_double_ptr_map.at(
id)->nu_pz = mydata.
data.Poissonpz;
564 if(mydata.
data.Shearzp!=0) {
565 shear_zp = mydata.
data.Shearzp;
567 shear_zp = mydata.
data.Youngz/(2*(1+mydata.
data.Poissonpz));
569 tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
570 tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
574 meshset,MBTET,trans_elastic.
setOfBlocks[
id].tEts,
true
581 trans_elastic.
setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(
id);
582 trans_elastic.
setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(
id);
586 if(trans_iso_blocks) {
597 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.
setOfBlocks.begin();
603 if(trans_iso_blocks) {
616 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.
setOfBlocks.begin();
620 "DISPLACEMENT",sit->second,trans_elastic.
commonData,2,
false,
false,
true
625 "DISPLACEMENT",sit->second,trans_elastic.
commonData
646 "DISPLACEMENT",sit->second,trans_elastic.
commonData,2,
true,
false,
true
651 "DISPLACEMENT",
"DISPLACEMENT",sit->second,trans_elastic.
commonData
660 "LAGRANGE_ELEM",
"DISPLACEMENT",
"LAGRANGE_MUL_DISP",
"MESH_NODE_POSITIONS"
669 "LAGRANGE_ELEM",
"DISPLACEMENT",
"LAGRANGE_MUL_TRAC",
"MESH_NODE_POSITIONS"
672 "LAGRANGE_ELEM_TRANS",
"DISPLACEMENT",
"LAGRANGE_MUL_RIGID_TRANS",
"MESH_NODE_POSITIONS"
675 "LAGRANGE_ELEM_ROT",
"DISPLACEMENT",
"LAGRANGE_MUL_RIGID_ROT",
"MESH_NODE_POSITIONS"
682 "LAGRANGE_ELEM",
"DISPLACEMENT",
"LAGRANGE_MUL_PERIODIC",
"MESH_NODE_POSITIONS",preriodic_prisms
685 "LAGRANGE_ELEM_TRANS",
"DISPLACEMENT",
"LAGRANGE_MUL_RIGID_TRANS",
"MESH_NODE_POSITIONS"
690 enum MINAMX { C0,MAXLAST };
692 ublas::vector<int> rowIndices;
695 rowIndices.resize(3*MAXLAST);
696 for(
int ii = 0;ii<6;ii++) {
697 rHs[ii].resize(3*MAXLAST);
701 MinMaxNodes minMaxNodes;
710 rval = moab.create_meshset(MESHSET_TRACK_OWNER,condensed_traction_element_meshset); CHKERRQ(
ierr);
713 if(it->getName().compare(0,12,
"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
719 nodes.merge(tris_nodes);
725 x.resize(nodes.size(),
false);
726 y.resize(nodes.size(),
false);
727 z.resize(nodes.size(),
false);
730 for(
int sx = -1; sx<=+1; sx+=2) {
731 for(
int sy = -1; sy<=+1; sy+=2) {
732 for(
int sz = -1; sz<=+1; sz+=2) {
733 if(bc_nb == MinMaxNodes::MAXLAST)
break;
735 dist_up_right.resize(x.size(),
false);
736 dist_up_right.clear();
737 for(
unsigned int nn = 0;nn<x.size();nn++) {
743 dist_up_right[nn] = sx*x[nn]+sy*y[nn]+sz*z[nn];
746 VectorDouble::iterator dist_it;
747 dist_it = max_element(dist_up_right.begin(),dist_up_right.end());
748 minMaxNodes.entMinMax[bc_nb++] = nodes[std::distance(dist_up_right.begin(),dist_it)];
765 if(trans_iso_blocks) {
767 "ELASTIC_MECHANICS",
"TRAN_ISOTROPIC_ELASTIC"
789 ierr = m_field.build_problems(); CHKERRQ(
ierr);
795 ierr = m_field.partition_problem(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
796 ierr = m_field.partition_finite_elements(
800 ierr = m_field.partition_ghost_dofs(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
805 ierr = m_field.VecCreateGhost(
"ELASTIC_MECHANICS",
ROW,&F[0]); CHKERRQ(
ierr);
806 for(
int ii = 1;ii<6;ii++) {
807 ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(
ierr);
809 ierr = m_field.VecCreateGhost(
"ELASTIC_MECHANICS",
COL,&
D); CHKERRQ(
ierr);
812 ierr = m_field.MatCreateMPIAIJWithArrays(
"ELASTIC_MECHANICS",&Aij); CHKERRQ(
ierr);
813 ierr = MatSetOption(Aij,MAT_STRUCTURALLY_SYMMETRIC,PETSC_TRUE); CHKERRQ(
ierr);
814 ierr = MatSetOption(Aij,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE); CHKERRQ(
ierr);
815 ierr = MatSetOption(Aij,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE); CHKERRQ(
ierr);
816 ierr = MatSetOption(Aij,MAT_USE_INODES,PETSC_TRUE); CHKERRQ(
ierr);
817 ierr = MatSetOption(Aij,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE); CHKERRQ(
ierr);
818 ierr = MatSetOption(Aij,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE); CHKERRQ(
ierr);
827 ierr = VecGhostUpdateBegin(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
828 ierr = VecGhostUpdateEnd(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
829 ierr = m_field.set_global_ghost_vector(
830 "ELASTIC_MECHANICS",
ROW,
D,INSERT_VALUES,SCATTER_REVERSE
832 for(
int ii = 0;ii<6;ii++) {
833 ierr = VecZeroEntries(F[ii]); CHKERRQ(
ierr);
834 ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
835 ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
837 ierr = MatZeroEntries(Aij); CHKERRQ(
ierr);
1053 int volume_vec_ghost[] = { 0 };
1054 ierr = VecCreateGhost(
1055 PETSC_COMM_WORLD,(!m_field.
get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
1057 ierr = VecZeroEntries(volume_vec); CHKERRQ(
ierr);
1072 lagrangian_element_disp.
setRVEBCsOperators(
"DISPLACEMENT",
"LAGRANGE_MUL_DISP",
"MESH_NODE_POSITIONS",Aij,F);
1077 lagrangian_element_trac.
setRVEBCsOperators(
"DISPLACEMENT",
"LAGRANGE_MUL_TRAC",
"MESH_NODE_POSITIONS",Aij,F);
1081 "DISPLACEMENT",
"LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_trac.
setOfRVEBC
1084 "ELASTIC_MECHANICS",
"LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.
getLoopFeRVEBCLhs()
1087 "DISPLACEMENT",
"LAGRANGE_MUL_RIGID_ROT",Aij,lagrangian_element_trac.
setOfRVEBC
1090 "ELASTIC_MECHANICS",
"LAGRANGE_ELEM_ROT",lagrangian_element_rigid_body_rot.
getLoopFeRVEBCLhs()
1094 lagrangian_element_periodic.
setRVEBCsOperators(
"DISPLACEMENT",
"LAGRANGE_MUL_PERIODIC",
"MESH_NODE_POSITIONS",Aij,F);
1096 "ELASTIC_MECHANICS",
"LAGRANGE_ELEM",lagrangian_element_periodic.
getLoopFeRVEBCLhs()
1099 "ELASTIC_MECHANICS",
"LAGRANGE_ELEM",lagrangian_element_periodic.
getLoopFeRVEBCRhs()
1102 "DISPLACEMENT",
"LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_periodic.
setOfRVEBC
1105 "ELASTIC_MECHANICS",
"LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.
getLoopFeRVEBCLhs()
1109 for(
int ii = 0;ii<6;ii++) {
1110 ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(
ierr);
1111 ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(
ierr);
1112 ierr = VecAssemblyBegin(F[ii]); CHKERRQ(
ierr);
1113 ierr = VecAssemblyEnd(F[ii]); CHKERRQ(
ierr);
1115 ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(
ierr);
1116 ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(
ierr);
1171 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(
ierr);
1172 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(
ierr);
1174 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(
ierr);
1175 ierr = PetscPrintf(PETSC_COMM_WORLD,
"RVE Volume %3.2g\n",rve_volume); CHKERRQ(
ierr);
1177 ierr = VecDestroy(&volume_vec);
1181 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(
ierr);
1182 ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(
ierr);
1183 ierr = KSPSetFromOptions(solver); CHKERRQ(
ierr);
1184 ierr = KSPSetUp(solver); CHKERRQ(
ierr);
1192 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
1195 PetscBool stress_by_boundary_integral = PETSC_FALSE;
1196 ierr = VecCreateGhost(
1197 PETSC_COMM_WORLD,(!m_field.
get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
1199 switch(choise_value) {
1202 "DISPLACEMENT",
"LAGRANGE_MUL_DISP",
"MESH_NODE_POSITIONS",stress_homo
1207 "DISPLACEMENT",
"LAGRANGE_MUL_TRAC",
"MESH_NODE_POSITIONS",stress_homo
1212 "DISPLACEMENT",
"LAGRANGE_MUL_PERIODIC",
"MESH_NODE_POSITIONS",stress_homo
1363 void setDoPreProcess() { doPreProcess =
true; }
1364 void unSetDoPreProcess() { doPreProcess =
false; }
1365 void setDoPostProcess() { doPostProcess =
true; }
1366 void unSetDoPostProcess() { doPostProcess =
false; }
1370 PetscErrorCode preProcess() {
1375 PetscFunctionReturn(0);
1377 PetscErrorCode postProcess() {
1382 PetscFunctionReturn(0);
1389 ierr = post_proc.addFieldValuesPostProc(
"DISPLACEMENT"); CHKERRQ(
ierr);
1390 ierr = post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT"); CHKERRQ(
ierr);
1391 ierr = post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS"); CHKERRQ(
ierr);
1392 ierr = post_proc.addFieldValuesGradientPostProc(
"MESH_NODE_POSITIONS"); CHKERRQ(
ierr);
1393 if(trans_iso_blocks) {
1394 ierr = post_proc.addFieldValuesGradientPostProc(
"POTENTIAL_FIELD"); CHKERRQ(
ierr);
1397 map<int,NonlinearElasticElement::BlockData>::iterator sit = iso_elastic.
setOfBlocks.begin();
1400 post_proc.getOpPtrVector().push_back(
1402 post_proc.postProcMesh,
1403 post_proc.mapGaussPts,
1406 post_proc.commonData,
1412 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.
setOfBlocks.begin();
1415 post_proc.getOpPtrVector().push_back(
1417 post_proc.postProcMesh,
1418 post_proc.mapGaussPts,
1421 post_proc.commonData
1427 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(
ierr);
1428 for(
int ii = 0;ii<6;ii++) {
1429 ierr = VecZeroEntries(
D); CHKERRQ(
ierr);
1430 ierr = KSPSolve(solver,F[ii],
D); CHKERRQ(
ierr);
1431 ierr = VecGhostUpdateBegin(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
1432 ierr = VecGhostUpdateEnd(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
1433 ierr = m_field.set_global_ghost_vector(
1434 "ELASTIC_MECHANICS",
ROW,
D,INSERT_VALUES,SCATTER_REVERSE
1436 post_proc.setDoPreProcess();
1437 post_proc.unSetDoPostProcess();
1439 "ELASTIC_MECHANICS",
"ELASTIC",post_proc
1441 post_proc.unSetDoPreProcess();
1442 post_proc.setDoPostProcess();
1446 "ELASTIC_MECHANICS",
"TRAN_ISOTROPIC_ELASTIC",post_proc
1450 sss <<
"mode_" <<
homo_bc_names[choise_value] <<
"_" << ii <<
".h5m";
1451 rval = post_proc.postProcMesh.write_file(
1452 sss.str().c_str(),
"MOAB",
"PARALLEL=WRITE_PART"
1455 ierr = VecZeroEntries(stress_homo); CHKERRQ(
ierr);
1495 PETSC_NULL,
"-my_rve_volume",&rve_volume,PETSC_NULL
1497 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(
ierr);
1498 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(
ierr);
1499 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(
ierr);
1500 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
1501 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
1502 for(
int jj=0; jj<6; jj++) {
1503 Dmat(jj,ii) = avec[jj];
1506 ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(
ierr);
1509 PETSC_COMM_WORLD,
"\nHomogenised Stiffens Matrix = \n\n"
1512 for(
int ii=0; ii<6; ii++) {
1515 "stress %d\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\n",
1516 ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
1523 char output_file_Dmat[255];
1528 if(pcomm->rank()==0){
1530 PetscViewer view_out;
1531 PetscViewerBinaryOpen(PETSC_COMM_SELF,output_file_Dmat,FILE_MODE_WRITE,&view_out);
1532 PetscViewerBinaryGetDescriptor(view_out,&fd);
1533 PetscBinaryWrite(fd,&Dmat(0,0),36,PETSC_DOUBLE,PETSC_FALSE);
1534 PetscViewerDestroy(&view_out);
1552 for(
int ii = 0;ii<6;ii++) {
1553 ierr = VecDestroy(&F[ii]); CHKERRQ(
ierr);
1556 ierr = MatDestroy(&Aij); CHKERRQ(
ierr);
1557 ierr = KSPDestroy(&solver); CHKERRQ(
ierr);
1558 ierr = VecDestroy(&stress_homo); CHKERRQ(
ierr);
1560 ierr = PetscTime(&v2);CHKERRQ(
ierr);
1561 ierr = PetscGetCPUTime(&t2);CHKERRQ(
ierr);
1563 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
Implementation of linear elastic material.
Operators and data structures for non-linear elastic analysis.
Post-process fields on refined mesh.
Post-processing stresses for non-linear analysis.
DEPRECATED typedef PostProcStress PostPorcStress
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
@ NOFIELD
scalar or vector of scalars describe (no true field)
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
#define CHKERRQ_MOAB(a)
check error code of MoAB function
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > >> > Face_CenPos_Handle_multiIndex
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
const double D
diffusivity
multi_index_container< Face_CenPos_Handle, indexed_by< ordered_non_unique< tag< xcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord > >, ordered_non_unique< tag< ycoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord > >, ordered_non_unique< tag< zcoord_tag >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > >, ordered_unique< tag< Tri_Hand_tag >, member< Face_CenPos_Handle, const EntityHandle,&Face_CenPos_Handle::Tri_Hand > >, ordered_unique< tag< Composite_xyzcoord >, composite_key< Face_CenPos_Handle, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::xcoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::ycoord >, member< Face_CenPos_Handle, double,&Face_CenPos_Handle::zcoord > > >> > Face_CenPos_Handle_multiIndex
int main(int argc, char *argv[])
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
const char * homo_bc_names[]
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
MyTriFE & getLoopFeRVEBCRhs()
MyTriFE & getLoopFeRVEBCStress()
MyTriFE & getLoopFeRVEBCLhs()
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec Stress_Homo)
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec)
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
MyPrismFE & getLoopFeRVEBCRhs()
MyPrismFE & getLoopFeRVEBCLhs()
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
MyPrismFE & getLoopFeRVEBCStress()
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
PetscErrorCode setRVEBCsRigidBodyRotOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
PetscErrorCode setRVEBCsRigidBodyTranOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
virtual int get_comm_size() const =0
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
virtual int get_comm_rank() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Transverse Isotropic material data structure.
Elastic material data structure.
Projection of edge entities with one mid-node on hierarchical basis.
Mat & snes_B
preconditioner of jacobian matrix
MoFEMErrorCode generateReferenceElementMesh()
definition of volume element
Operator performs automatic differentiation.
structure grouping operators and data used for calculation of nonlinear elastic element
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MyVolumeFE & getLoopFeLhs()
get lhs volume element
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop