v0.13.1
Classes | Typedefs | Enumerations | Functions | Variables
rve_mechanical.cpp File Reference

Calculates stiffness matrix for elastic RVE. More...

#include <MoFEM.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <BCs_RVELagrange_Disp.hpp>
#include <BCs_RVELagrange_Trac.hpp>
#include <BCs_RVELagrange_Trac_Rigid_Rot.hpp>
#include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
#include <BCs_RVELagrange_Periodic.hpp>
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <Hooke.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <SmallTransverselyIsotropic.hpp>
#include <VolumeCalculation.hpp>
#include <boost/ptr_container/ptr_map.hpp>
#include <PostProcOnRefMesh.hpp>
#include <PostProcStresses.hpp>
#include <moab/AdaptiveKDTree.hpp>
#include <algorithm>

Go to the source code of this file.

Classes

struct  Face_CenPos_Handle
 
struct  xcoord_tag
 
struct  ycoord_tag
 
struct  zcoord_tag
 
struct  Tri_Hand_tag
 
struct  Composite_xyzcoord
 

Typedefs

typedef 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
 

Enumerations

enum  HomoBCTypes {
  HOMOBCDISP , HOMOBCPERIODIC , HOMOBCTRAC , NITSCHE ,
  NOHOMOBC
}
 

Functions

void tetcircumcenter_tp (double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
 
void tricircumcenter3d_tp (double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
 
int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
const char * homo_bc_names []
 

Detailed Description

Calculates stiffness matrix for elastic RVE.

Three types of boundary conditions are implemented, i.e. HOMOBCDISP, HOMOBCPERIODIC, HOMOBCTRAC, NITSCHE.

NITSHCE method allow to apply periodic boundary conditions to arbitrary convex shape RVE.

Definition in file rve_mechanical.cpp.

Typedef Documentation

◆ Face_CenPos_Handle_multiIndex

typedef 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

Definition at line 134 of file rve_mechanical.cpp.

Enumeration Type Documentation

◆ HomoBCTypes

enum HomoBCTypes
Enumerator
HOMOBCDISP 
HOMOBCPERIODIC 
HOMOBCTRAC 
NITSCHE 
NOHOMOBC 

Definition at line 79 of file rve_mechanical.cpp.

79 {
80 HOMOBCDISP,
81 HOMOBCPERIODIC,
82 HOMOBCTRAC,
83 NITSCHE,
84 NOHOMOBC
85};
NOHOMOBC
@ NOHOMOBC
Definition: rve_mechanical.cpp:84
HOMOBCPERIODIC
@ HOMOBCPERIODIC
Definition: rve_mechanical.cpp:81
NITSCHE
@ NITSCHE
Definition: rve_mechanical.cpp:83
HOMOBCDISP
@ HOMOBCDISP
Definition: rve_mechanical.cpp:80
HOMOBCTRAC
@ HOMOBCTRAC
Definition: rve_mechanical.cpp:82

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 136 of file rve_mechanical.cpp.

136 {
137
138 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
139
140 try {
141
142 moab::Core mb_instance;
143 moab::Interface& moab = mb_instance;
144 int rank;
145 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
146
147 //Reade parameters from line command
148 PetscBool flg = PETSC_TRUE;
149 char mesh_file_name[255];
150 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
151 if(flg != PETSC_TRUE) {
152 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_FOUND,"*** ERROR -my_file (MESH FILE NEEDED)");
153 }
154
155 PetscInt order;
156 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
157 if(flg != PETSC_TRUE) {
158 order = 1;
159 }
160
161 PetscInt choise_value = NOHOMOBC;
162 ierr = PetscOptionsGetEList(
163 NULL,"-my_bc_type",homo_bc_names,NOHOMOBC,&choise_value,&flg
164 ); CHKERRQ(ierr);
165 if(flg != PETSC_TRUE) {
166 SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"*** Boundary conditions not set");
167 }
168
169 //Read mesh to MOAB
170 const char *option;
171 option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
172 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
173 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
174 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
175
176 //We need that for code profiling
177 PetscLogDouble t1,t2;
178 PetscLogDouble v1,v2;
179 ierr = PetscTime(&v1); CHKERRQ(ierr);
180 ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
181
182 //Create MoFEM (Joseph) database
183 MoFEM::Core core(moab);
184 MoFEM::Interface& m_field = core;
185
186 vector<BitRefLevel> bit_levels;
187 {
188 Tag th_meshset_info;
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
192 );
193 int meshset_data[2];
194 EntityHandle root = moab.get_root_set();
195 rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
196 if(meshset_data[0]==0) {
197 meshset_data[0] = 1;
198 rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
199
200 }
201 bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
202 }
203 BitRefLevel problem_bit_level = bit_levels.back();
204 ierr = m_field.seed_ref_level_3D(0,problem_bit_level); CHKERRQ(ierr);
205
206 // const clock_t begin_time = clock();
207 ierr = m_field.build_fields(); CHKERRQ(ierr);
208 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
209
210 Range preriodic_prisms;
211 if(choise_value == HOMOBCPERIODIC) {
212 //FIXME: Naming convention is not consistent in this section of code
213 //=======================================================================================================
214 //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
215 //=======================================================================================================
216 //Populating the Multi-index container with -ve triangles
217 Range SurTrisNeg;
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);
220 Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
221 double TriCen[3], coords_Tri[9];
222
223
224 double roundfact=1e3;
225
226 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
227 const EntityHandle* conn_face; int num_nodes_Tri;
228
229 //get nodes attached to the face
230 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
231 //get nodal coordinates
232 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
233
234 //Find triangle centriod
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;
238
239 //round values to roundfact disimal places
240 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
241 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
242 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
243
244 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
245 //fill the multi-index container with centriod coordinates and triangle handles
246 Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
247 // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
248 // cout<<endl<<endl;
249 }
250
251 //Populating the Multi-index container with +ve triangles
252 Range SurTrisPos;
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++) {
256 const EntityHandle* conn_face; int num_nodes_Tri;
257
258 //get nodes attached to the face
259 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
260 //get nodal coordinates
261 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
262
263 //Find triangle centriod
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;
267
268 //round values to roundfact disimal places
269 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
270 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
271 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
272 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
273
274 //fill the multi-index container with centriod coordinates and triangle handles
275 Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
276 }
277
278 //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
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;
286
287 //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
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;
294
295 /*double LxRVE, LyRVE, LzRVE;
296 LxRVE=XcoordMax-XcoordMin;
297 LyRVE=YcoordMax-YcoordMin;
298 LzRVE=ZcoordMax-ZcoordMin;*/
299
300 //Creating Prisims between triangles on -ve and +ve faces
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;
306 //int count=0;
307
308 //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
309 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
310
311 Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
312 // cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
313
314 //corresponding +ve triangle
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));
319 // cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
320 // cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
321
322
323 //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
324 EntityHandle PrismNodes[6];
325 vector<EntityHandle> TriNodesNeg, TriNodesPos;
326 double CoordNodeNeg[9], CoordNodePos[9];
327 rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
328 rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
329 rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
330 rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
331 for(int ii=0; ii<3; ii++){
332 PrismNodes[ii]=TriNodesNeg[ii];
333 }
334 // for(int ii=0; ii<3; ii++){
335 // cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
336 // }
337 // for(int ii=0; ii<3; ii++){
338 // cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
339 // }
340
341 //Match exact nodes to each other to avoide the problem of twisted prisms
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++){
348 //Round nodal coordinates to 3 dicimal places only for comparison //round values to 3 disimal places
349 XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
350 YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
351 ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
352
353
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;
358
359 if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
360 PrismNodes[3+ii]=TriNodesPos[jj];
361 break;
362 }
363 }
364 }
365 //prism nodes and their coordinates
366 double CoordNodesPrisms[18];
367 rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
368 // for(int ii=0; ii<6; ii++){
369 // cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
370 // }
371 // cout<<endl<<endl;
372 //insertion of individula prism element and its addition to range preriodic_prisms
373 EntityHandle PeriodicPrism;
374 rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
375 preriodic_prisms.insert(PeriodicPrism);
376
377 }
378
379 //insertion of individual prism element and its addition to range preriodic_prisms
380 ierr = m_field.seed_ref_level(preriodic_prisms,problem_bit_level); CHKERRQ(ierr);
381
382 }
383
384 EntityHandle out_meshset;
385 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
386 // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
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);
392
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
396 ); CHKERRQ(ierr);
397 //to create meshset from range
398 EntityHandle meshset_prims_on_problem_bit_level;
399 rval = moab.create_meshset(MESHSET_SET,meshset_prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
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);
402
403 //Fields
404 int field_rank=3;
405 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
406 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
407 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
408 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
409
410 //add entitities (by tets) to the field
411 if(choise_value == HOMOBCDISP) {
412 ierr = m_field.add_field("LAGRANGE_MUL_DISP",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
413 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
414 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
415 Range tris;
416 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
417 ierr = m_field.add_ents_to_field_by_type(tris,MBTRI,"LAGRANGE_MUL_DISP"); CHKERRQ(ierr);
418 }
419 }
420 }
421
422 if(choise_value == HOMOBCPERIODIC) {
423 ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
424 //Control 3 rigid body translations in x, y and z axis
425 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr);
426 // Setting up dummy no-field vertex
427 EntityHandle no_field_vertex;
428 {
429 const double coords[] = {0,0,0};
430 rval = m_field.get_moab().create_vertex(coords,no_field_vertex); CHKERRQ_MOAB(rval);
431 Range range_no_field_vertex;
432 range_no_field_vertex.insert(no_field_vertex);
433 ierr = m_field.seed_ref_level(range_no_field_vertex,BitRefLevel().set()); CHKERRQ(ierr);
434 EntityHandle meshset;
435 meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_TRANS");
436 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
437 }
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);
441 ierr = m_field.add_ents_to_field_by_type(surface_negative,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
442 }
443
444 if(choise_value == HOMOBCTRAC) {
445 ierr = m_field.add_field("LAGRANGE_MUL_TRAC",NOFIELD,NOBASE,6); CHKERRQ(ierr);
446 //Control 3 rigid body translations in x, y and z axis
447 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr);
448 //Controla 3 rigid body rotations about x, y and z axis
449 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_ROT",NOFIELD,NOBASE,3); CHKERRQ(ierr);
450 EntityHandle no_field_vertex;
451 {
452 const double coords[] = {0,0,0};
453 rval = m_field.get_moab().create_vertex(coords,no_field_vertex); CHKERRQ_MOAB(rval);
454 Range range_no_field_vertex;
455 range_no_field_vertex.insert(no_field_vertex);
456 ierr = m_field.seed_ref_level(range_no_field_vertex,BitRefLevel().set()); CHKERRQ(ierr);
457 EntityHandle meshset;
458 meshset = m_field.get_field_meshset("LAGRANGE_MUL_TRAC");
459 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
460 meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_TRANS");
461 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
462 meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_ROT");
463 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
464 }
465 }
466
467 //set app. order
468 //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
469 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
470 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
471 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
472 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
473
474 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
475 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
476 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
477 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
478
479 PetscBool fo_boundary = PETSC_FALSE;
480 ierr = PetscOptionsGetBool(PETSC_NULL,"-my_fo_boundary",&fo_boundary,PETSC_NULL); CHKERRQ(ierr);
481 if(fo_boundary) {
482 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
483 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
484 Range tris;
485 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
486 Range tris_edges;
487 rval = moab.get_adjacencies(tris,1,false,tris_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
488 ierr = m_field.set_field_order(tris,"DISPLACEMENT",1); CHKERRQ(ierr);
489 ierr = m_field.set_field_order(tris_edges,"DISPLACEMENT",1); CHKERRQ(ierr);
490 }
491 }
492 }
493
494 if(choise_value == HOMOBCDISP) {
495 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
496 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
497 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_DISP",1); CHKERRQ(ierr);
498 }
499
500 if(choise_value == HOMOBCPERIODIC) {
501 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
502 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
503 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
504 }
505
506 //build field
507 ierr = m_field.build_fields(); CHKERRQ(ierr);
508
509 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
510 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
511
512 boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>());
513 boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>());
514
515 NonlinearElasticElement iso_elastic(m_field,1);
516 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|MAT_ELASTICSET,it)) {
517 if(it->getName() != "MAT_ELASTIC_1") continue;
518 Mat_Elastic mydata;
519 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
520 int id = it->getMeshsetId();
521 EntityHandle meshset = it->getMeshset();
522 rval = m_field.get_moab().get_entities_by_type(
523 meshset,MBTET,iso_elastic.setOfBlocks[id].tEts,true
524 ); CHKERRQ_MOAB(rval);
525 iso_elastic.setOfBlocks[id].iD = id;
526 iso_elastic.setOfBlocks[id].E = mydata.data.Young;
527 iso_elastic.setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
528 iso_elastic.setOfBlocks[id].materialDoublePtr = hooke_double;
529 iso_elastic.setOfBlocks[id].materialAdoublePtr = hooke_adouble;
530 ierr = m_field.seed_finite_elements(iso_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
531 }
532 ierr = iso_elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
533 ierr = iso_elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
534 if(m_field.check_field("POTENTIAL_FIELD")) {
535 ierr = m_field.modify_finite_element_add_field_data("ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
536 }
537
538 NonlinearElasticElement trans_elastic(m_field,2);
539 trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
540 trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
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;
544 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
545 //Get block name
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();
550 Mat_Elastic_TransIso mydata;
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>();
554 //nu_p, nu_pz, E_p, E_z, G_zp
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;
563 double shear_zp;
564 if(mydata.data.Shearzp!=0) {
565 shear_zp = mydata.data.Shearzp;
566 } else {
567 shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
568 }
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;
571 //get tets from block where material is defined
572 EntityHandle meshset = it->getMeshset();
573 rval = m_field.get_moab().get_entities_by_type(
574 meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
575 ); CHKERRQ_MOAB(rval);
576 //adding material to nonlinear class
577 trans_elastic.setOfBlocks[id].iD = id;
578 //note that material parameters are defined internally in material model
579 trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
580 trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
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);
583 ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
584 }
585 }
586 if(trans_iso_blocks) {
587 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
588 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
589 ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
590 ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
591 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
592 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
593 if(m_field.check_field("MESH_NODE_POSITIONS")) {
594 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
595 }
596 for(
597 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
598 sit!=trans_elastic.setOfBlocks.end();sit++
599 ) {
600 ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
601 }
602 }
603 if(trans_iso_blocks) {
604 //Rhs
605 trans_elastic.feRhs.getOpPtrVector().push_back(
606 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData)
607 );
608 trans_elastic.feRhs.getOpPtrVector().push_back(
609 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData)
610 );
611 if(m_field.check_field("MESH_NODE_POSITIONS")) {
612 trans_elastic.feRhs.getOpPtrVector().push_back(
613 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData)
614 );
615 }
616 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
617 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
618 trans_elastic.feRhs.getOpPtrVector().push_back(
619 new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
620 "DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true
621 )
622 );
623 trans_elastic.feRhs.getOpPtrVector().push_back(
624 new NonlinearElasticElement::OpRhsPiolaKirchhoff(
625 "DISPLACEMENT",sit->second,trans_elastic.commonData
626 )
627 );
628 }
629
630 //Lhs
631 trans_elastic.feLhs.getOpPtrVector().push_back(
632 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData)
633 );
634 trans_elastic.feLhs.getOpPtrVector().push_back(
635 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData)
636 );
637 if(m_field.check_field("MESH_NODE_POSITIONS")) {
638 trans_elastic.feLhs.getOpPtrVector().push_back(
639 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData)
640 );
641 }
642 sit = trans_elastic.setOfBlocks.begin();
643 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
644 trans_elastic.feLhs.getOpPtrVector().push_back(
645 new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
646 "DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true
647 )
648 );
649 trans_elastic.feLhs.getOpPtrVector().push_back(
650 new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(
651 "DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData
652 )
653 );
654 }
655 }
656
657 BCs_RVELagrange_Disp lagrangian_element_disp(m_field);
658 if(choise_value == HOMOBCDISP) {
659 lagrangian_element_disp.addLagrangiangElement(
660 "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS"
661 );
662 }
663
664 BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
665 BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
666 BCs_RVELagrange_Trac_Rigid_Rot lagrangian_element_rigid_body_rot(m_field);
667 if(choise_value == HOMOBCTRAC) {
668 lagrangian_element_trac.addLagrangiangElement(
669 "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS"
670 );
671 lagrangian_element_trac.addLagrangiangElement(
672 "LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS"
673 );
674 lagrangian_element_trac.addLagrangiangElement(
675 "LAGRANGE_ELEM_ROT","DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT","MESH_NODE_POSITIONS"
676 );
677 }
678
679 BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
680 if(choise_value == HOMOBCPERIODIC) {
681 lagrangian_element_periodic.addLagrangiangElement(
682 "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",preriodic_prisms
683 );
684 lagrangian_element_trac.addLagrangiangElement(
685 "LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS"
686 );
687 }
688
689 struct MinMaxNodes {
690 enum MINAMX { C0,MAXLAST };
691 EntityHandle entMinMax[MAXLAST];
692 ublas::vector<int> rowIndices;
693 VectorDouble rHs[6];
694 MinMaxNodes() {
695 rowIndices.resize(3*MAXLAST);
696 for(int ii = 0;ii<6;ii++) {
697 rHs[ii].resize(3*MAXLAST);
698 }
699 }
700 };
701 MinMaxNodes minMaxNodes;
702
703 if(choise_value == NITSCHE) { // Condensed traction/periodc BC
704 ierr = m_field.add_finite_element("SURFACE_ELEMENTS"); CHKERRQ(ierr);
705 ierr = m_field.modify_finite_element_add_field_row("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
706 ierr = m_field.modify_finite_element_add_field_col("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
707 ierr = m_field.modify_finite_element_add_field_data("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
708 ierr = m_field.modify_finite_element_add_field_data("SURFACE_ELEMENTS","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
709 EntityHandle condensed_traction_element_meshset;
710 rval = moab.create_meshset(MESHSET_TRACK_OWNER,condensed_traction_element_meshset); CHKERRQ(ierr);
711 Range nodes;
712 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
713 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
714 Range tris;
715 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
716 ierr = m_field.add_ents_to_finite_element_by_type(tris,MBTRI,"SURFACE_ELEMENTS"); CHKERRQ(ierr);
717 Range tris_nodes;
718 rval = moab.get_connectivity(tris,nodes,true); CHKERRQ_MOAB(rval);
719 nodes.merge(tris_nodes);
720 }
721 }
722
723 {
724 VectorDouble x,y,z;
725 x.resize(nodes.size(),false);
726 y.resize(nodes.size(),false);
727 z.resize(nodes.size(),false);
728 rval = moab.get_coords(nodes,&x[0],&y[0],&z[0]); CHKERRQ_MOAB(rval);
729 int bc_nb = 0;
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;
734 VectorDouble dist_up_right;
735 dist_up_right.resize(x.size(),false);
736 dist_up_right.clear();
737 for(unsigned int nn = 0;nn<x.size();nn++) {
738 if(
739 ((sx*x[nn])>0)&&
740 ((sy*y[nn])>0)&&
741 ((sz*z[nn])>0)
742 ) {
743 dist_up_right[nn] = sx*x[nn]+sy*y[nn]+sz*z[nn];
744 }
745 }
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)];
749 }
750 }
751 }
752 }
753
754 }
755
756 //build finite elements
757 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
758 //build adjacencies
759 ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
760
761 //define problems
762 ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
763 //set finite elements for problem
764 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
765 if(trans_iso_blocks) {
766 ierr = m_field.modify_problem_add_finite_element(
767 "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC"
768 ); CHKERRQ(ierr);
769 }
770 if(choise_value == HOMOBCDISP) {
771 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
772 }
773 if(choise_value == HOMOBCTRAC) {
774 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
775 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
776 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT"); CHKERRQ(ierr);
777 }
778 if(choise_value == HOMOBCPERIODIC) {
779 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
780 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
781 }
782 if(choise_value == NITSCHE) {
783 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","SURFACE_ELEMENTS"); CHKERRQ(ierr);
784 }
785
786 //set refinement level for problem
787 ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",problem_bit_level); CHKERRQ(ierr);
788 //build problem
789 ierr = m_field.build_problems(); CHKERRQ(ierr);
790
791 /****/
792 //mesh partitioning
793
794 //partition
795 ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
796 ierr = m_field.partition_finite_elements(
797 "ELASTIC_MECHANICS",false,0,m_field.get_comm_size() // build elements on all procs
798 ); CHKERRQ(ierr);
799 //what are ghost nodes, see Petsc Manual
800 ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
801
802 //create matrices
803 Vec D;
804 vector<Vec> F(6);
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);
808 }
809 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
810
811 Mat Aij;
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);
819
820 /*{
821 ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
822 std::string wait;
823 std::cin >> wait;
824 }*/
825
826 ierr = VecZeroEntries(D); 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
831 ); CHKERRQ(ierr);
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);
836 }
837 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
838
839 // NitscheMethod::BlockData nitsche_block_data;
840 // NitscheMethod::CommonData nitsche_common_data;
841 // PeriodicNitscheConstrains::CommonData periodic_common_data;
842 // PeriodicNitscheConstrains::MyNitscheVolume nitsche_element_iso(
843 // m_field,nitsche_block_data,nitsche_common_data,periodic_common_data
844 // );
845 // PeriodicNitscheConstrains::MyNitscheVolume nitsche_element_trans(
846 // m_field,nitsche_block_data,nitsche_common_data,periodic_common_data
847 // );
848
849 // if(choise_value == NITSCHE) {
850
851 // nitsche_block_data.faceElemName = "SURFACE_ELEMENTS";
852 // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
853 // if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
854 // Range tris;
855 // rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
856 // nitsche_block_data.fAces.merge(tris);
857 // }
858 // }
859
860 // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
861 // if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
862 // Range tris;
863 // rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
864 // periodic_common_data.skinFaces.merge(tris);
865 // }
866 // }
867
868 // nitsche_block_data.gamma = 1e-7;
869 // nitsche_block_data.phi = -1;
870 // periodic_common_data.eps = 0;
871 // ierr = PetscOptionsGetReal(
872 // PETSC_NULL,"-my_gamma",&nitsche_block_data.gamma,PETSC_NULL
873 // ); CHKERRQ(ierr);
874 // ierr = PetscOptionsGetReal(
875 // PETSC_NULL,"-my_phi",&nitsche_block_data.phi,PETSC_NULL
876 // ); CHKERRQ(ierr);
877 // ierr = PetscOptionsGetReal(
878 // PETSC_NULL,"-my_eps",&periodic_common_data.eps,PETSC_NULL
879 // ); CHKERRQ(ierr);
880 // ierr = PetscPrintf(
881 // PETSC_COMM_WORLD,
882 // "Nitsche method gamma = %4.2e phi = %2.1f eps = %4.2e\n",
883 // nitsche_block_data.gamma,nitsche_block_data.phi,periodic_common_data.eps
884 // ); CHKERRQ(ierr);
885
886 // for(
887 // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
888 // mit!=iso_elastic.setOfBlocks.end();
889 // mit++
890 // ) {
891 // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
892 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
893 // nitsche_element_iso.snes_B = Aij;
894
895 // nitsche_element_iso.getOpPtrVector().push_back(
896 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
897 // "DISPLACEMENT",elastic_common_data
898 // )
899 // );
900 // nitsche_element_iso.getOpPtrVector().push_back(
901 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
902 // "MESH_NODE_POSITIONS",elastic_common_data
903 // )
904 // );
905 // nitsche_element_iso.getOpPtrVector().push_back(
906 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
907 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,true,false,true
908 // )
909 // );
910 // nitsche_element_iso.getOpPtrVector().push_back(
911 // new PeriodicNitscheConstrains::OpLhsPeriodicNormal(
912 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
913 // elastic_block_data,elastic_common_data,
914 // periodic_common_data
915 // )
916 // );
917 // nitsche_element_iso.getOpPtrVector().push_back(
918 // new PeriodicNitscheConstrains::OpRhsPeriodicNormal(
919 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
920 // elastic_block_data,elastic_common_data,
921 // periodic_common_data,
922 // F
923 // )
924 // );
925 // nitsche_element_iso.getOpPtrVector().push_back(
926 // new NitscheMethod::OpLhsNormal(
927 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
928 // elastic_block_data,elastic_common_data,true
929 // )
930 // );
931
932 // // this is to get data on opposite element
933 // nitsche_element_iso.periodicVolume.getOpPtrVector().clear();
934 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
935 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
936 // "DISPLACEMENT",elastic_common_data
937 // )
938 // );
939 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
940 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
941 // "MESH_NODE_POSITIONS",elastic_common_data
942 // )
943 // );
944 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
945 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
946 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,true,false,true
947 // )
948 // );
949 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
950 // new PeriodicNitscheConstrains::OpGetVolumeData(
951 // elastic_common_data,
952 // periodic_common_data
953 // )
954 // );
955 // periodic_common_data.volumeElemName = "ELASTIC";
956
957 // nitsche_element_iso.addToRule = 1;
958 // ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",nitsche_element_iso); CHKERRQ(ierr);
959 // }
960
961 // for(
962 // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
963 // mit!=trans_elastic.setOfBlocks.end();
964 // mit++
965 // ) {
966 // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
967 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
968 // nitsche_element_trans.snes_B = Aij;
969
970 // nitsche_element_trans.getOpPtrVector().push_back(
971 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
972 // "DISPLACEMENT",elastic_common_data
973 // )
974 // );
975 // nitsche_element_trans.getOpPtrVector().push_back(
976 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
977 // "MESH_NODE_POSITIONS",elastic_common_data
978 // )
979 // );
980 // nitsche_element_trans.getOpPtrVector().push_back(
981 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
982 // "POTENTIAL_FIELD",elastic_common_data
983 // )
984 // );
985 // nitsche_element_trans.getOpPtrVector().push_back(
986 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
987 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,true,false,true
988 // )
989 // );
990 // nitsche_element_trans.getOpPtrVector().push_back(
991 // new PeriodicNitscheConstrains::OpLhsPeriodicNormal(
992 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
993 // elastic_block_data,elastic_common_data,
994 // periodic_common_data
995 // )
996 // );
997 // nitsche_element_trans.getOpPtrVector().push_back(
998 // new PeriodicNitscheConstrains::OpRhsPeriodicNormal(
999 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1000 // elastic_block_data,elastic_common_data,
1001 // periodic_common_data,
1002 // F
1003 // )
1004 // );
1005 // nitsche_element_trans.getOpPtrVector().push_back(
1006 // new NitscheMethod::OpLhsNormal(
1007 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1008 // elastic_block_data,elastic_common_data,true
1009 // )
1010 // );
1011
1012 // // this is to get data on opposite element
1013 // nitsche_element_trans.periodicVolume.getOpPtrVector().clear();
1014 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1015 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1016 // "DISPLACEMENT",elastic_common_data
1017 // )
1018 // );
1019 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1020 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1021 // "POTENTIAL_FIELD",elastic_common_data
1022 // )
1023 // );
1024 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1025 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1026 // "MESH_NODE_POSITIONS",elastic_common_data
1027 // )
1028 // );
1029 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1030 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1031 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,true,false,true
1032 // )
1033 // );
1034 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1035 // new PeriodicNitscheConstrains::OpGetVolumeData(
1036 // elastic_common_data,
1037 // periodic_common_data
1038 // )
1039 // );
1040 // periodic_common_data.volumeElemName = "TRAN_ISOTROPIC_ELASTIC";
1041
1042 // nitsche_element_trans.addToRule = 1;
1043 // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1044 // ierr = m_field.loop_finite_elements(
1045 // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",nitsche_element_trans
1046 // ); CHKERRQ(ierr);
1047 // }
1048
1049
1050 // }
1051
1052 Vec volume_vec;
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
1056 ); CHKERRQ(ierr);
1057 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
1058
1059 iso_elastic.getLoopFeLhs().getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1060 trans_elastic.getLoopFeLhs().getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1061
1062 //iso_elastic element matrix
1063 iso_elastic.getLoopFeLhs().snes_x = D;
1064 iso_elastic.getLoopFeLhs().snes_B = Aij;
1065 trans_elastic.getLoopFeLhs().snes_x = D;
1066 trans_elastic.getLoopFeLhs().snes_B = Aij;
1067 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",iso_elastic.getLoopFeLhs()); CHKERRQ(ierr);
1068 if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1069 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",trans_elastic.getLoopFeLhs()); CHKERRQ(ierr);
1070
1071 if(choise_value == HOMOBCDISP) {
1072 lagrangian_element_disp.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",Aij,F);
1073 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
1074 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
1075 }
1076 if(choise_value == HOMOBCTRAC) {
1077 lagrangian_element_trac.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",Aij,F);
1078 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
1079 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
1080 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators(
1081 "DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_trac.setOfRVEBC
1082 );
1083 ierr = m_field.loop_finite_elements(
1084 "ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()
1085 ); CHKERRQ(ierr);
1086 lagrangian_element_rigid_body_rot.setRVEBCsRigidBodyRotOperators(
1087 "DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT",Aij,lagrangian_element_trac.setOfRVEBC
1088 );
1089 ierr = m_field.loop_finite_elements(
1090 "ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT",lagrangian_element_rigid_body_rot.getLoopFeRVEBCLhs()
1091 ); CHKERRQ(ierr);
1092 }
1093 if(choise_value == HOMOBCPERIODIC) {
1094 lagrangian_element_periodic.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,F);
1095 ierr = m_field.loop_finite_elements(
1096 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCLhs()
1097 ); CHKERRQ(ierr);
1098 ierr = m_field.loop_finite_elements(
1099 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCRhs()
1100 ); CHKERRQ(ierr);
1101 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators(
1102 "DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_periodic.setOfRVEBC
1103 );
1104 ierr = m_field.loop_finite_elements(
1105 "ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()
1106 ); CHKERRQ(ierr);
1107 }
1108
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);
1114 }
1115 ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1116 ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1117
1118 {
1119 //ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
1120 //std::string wait;
1121 //std::cin >> wait;
1122 }
1123
1124 // if(choise_value == NITSCHE) {
1125 // if(periodic_common_data.eps==0) {
1126 // const Problem *problem_ptr;
1127 // ierr = m_field.get_problem("ELASTIC_MECHANICS",&problem_ptr); CHKERRQ(ierr);
1128 // for(int nn = 0;nn!=MinMaxNodes::MAXLAST;nn++) {
1129 // for(
1130 // _IT_NUMEREDDOF_ROW_BY_ENT_FOR_LOOP_(
1131 // problem_ptr,minMaxNodes.entMinMax[nn],dof
1132 // )
1133 // ) {
1134 // minMaxNodes.rowIndices[3*nn+dof->get()->getDofCoeffIdx()]
1135 // = dof->get()->getPetscGlobalDofIdx();
1136 // }
1137 // }
1138 // VectorDouble coords;
1139 // int nb_bcs = 3*MinMaxNodes::MAXLAST;
1140 // coords.resize(nb_bcs,false);
1141 // rval = moab.get_coords(&minMaxNodes.entMinMax[0],MinMaxNodes::MAXLAST,&coords[0]);
1142 // VectorDouble strain;
1143 // strain.resize(6,false);
1144 // MatrixDouble mat_d;
1145 // for(int ii = 0;ii<6;ii++) {
1146 // minMaxNodes.rHs[ii].clear();
1147 // strain.clear();
1148 // strain[ii] = 1;
1149 // for(int nn = 0;nn<MinMaxNodes::MAXLAST;nn++) {
1150 // PeriodicNitscheConstrains::OpRhsPeriodicNormal::calcualteDMatrix(
1151 // VectorAdaptor(3,ublas::shallow_array_adaptor<double>(3,&coords[3*nn])),
1152 // mat_d
1153 // );
1154 // VectorAdaptor rhs(3,ublas::shallow_array_adaptor<double>(3,&minMaxNodes.rHs[ii][3*nn]));
1155 // noalias(rhs) = prod(mat_d,strain);
1156 // }
1157 // }
1158 // for(int ii = 0;ii<6;ii++) {
1159 // ierr = VecSetValues(
1160 // F[ii],nb_bcs,&minMaxNodes.rowIndices[0],&minMaxNodes.rHs[ii][0],INSERT_VALUES
1161 // ); CHKERRQ(ierr);
1162 // ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
1163 // ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
1164 // }
1165 // ierr = MatZeroRows(
1166 // Aij,nb_bcs,&minMaxNodes.rowIndices[0],1,PETSC_NULL,PETSC_NULL
1167 // ); CHKERRQ(ierr);
1168 // }
1169 // }
1170
1171 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
1172 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
1173 double rve_volume;
1174 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
1175 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
1176
1177 ierr = VecDestroy(&volume_vec);
1178
1179 // Solver
1180 KSP solver;
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);
1185
1186 MatrixDouble Dmat;
1187 Dmat.resize(6,6);
1188 Dmat.clear();
1189
1190 //create a vector for 6 components of homogenized stress
1191 Vec stress_homo;
1192 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
1193 NonlinearElasticElement::MyVolumeFE ave_stress_iso(m_field);
1194 NonlinearElasticElement::MyVolumeFE ave_stress_trans(m_field);
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
1198 ); CHKERRQ(ierr);
1199 switch(choise_value) {
1200 case HOMOBCDISP:
1201 lagrangian_element_disp.setRVEBCsHomoStressOperators(
1202 "DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",stress_homo
1203 );
1204 break;
1205 case HOMOBCTRAC:
1206 lagrangian_element_trac.setRVEBCsHomoStressOperators(
1207 "DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",stress_homo
1208 );
1209 break;
1210 case HOMOBCPERIODIC:
1211 lagrangian_element_periodic.setRVEBCsHomoStressOperators(
1212 "DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo
1213 );
1214 break;
1215 case NITSCHE:
1216 // if(stress_by_boundary_integral) {
1217 // nitsche_element_iso.getOpPtrVector().clear();
1218 // nitsche_element_trans.getOpPtrVector().clear();
1219 // nitsche_element_iso.periodicVolume.getOpPtrVector().clear();
1220 // nitsche_element_trans.periodicVolume.getOpPtrVector().clear();
1221 // for(
1222 // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
1223 // mit!=iso_elastic.setOfBlocks.end();
1224 // mit++
1225 // ) {
1226 // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
1227 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1228 // nitsche_element_iso.getOpPtrVector().push_back(
1229 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1230 // "DISPLACEMENT",elastic_common_data
1231 // )
1232 // );
1233 // nitsche_element_iso.getOpPtrVector().push_back(
1234 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1235 // "MESH_NODE_POSITIONS",elastic_common_data
1236 // )
1237 // );
1238 // nitsche_element_iso.getOpPtrVector().push_back(
1239 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1240 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,false,false,true
1241 // )
1242 // );
1243 // nitsche_element_iso.getOpPtrVector().push_back(
1244 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral(
1245 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1246 // elastic_block_data,elastic_common_data,stress_homo
1247 // )
1248 // );
1249 // }
1250 // for(
1251 // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
1252 // mit!=trans_elastic.setOfBlocks.end();
1253 // mit++
1254 // ) {
1255 // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
1256 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1257 // nitsche_element_trans.getOpPtrVector().push_back(
1258 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1259 // "DISPLACEMENT",elastic_common_data
1260 // )
1261 // );
1262 // nitsche_element_trans.getOpPtrVector().push_back(
1263 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1264 // "MESH_NODE_POSITIONS",elastic_common_data
1265 // )
1266 // );
1267 // nitsche_element_trans.getOpPtrVector().push_back(
1268 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1269 // "POTENTIAL_FIELD",elastic_common_data
1270 // )
1271 // );
1272 // nitsche_element_trans.getOpPtrVector().push_back(
1273 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1274 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,false,false,true
1275 // )
1276 // );
1277 // nitsche_element_trans.getOpPtrVector().push_back(
1278 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral(
1279 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1280 // elastic_block_data,elastic_common_data,stress_homo
1281 // )
1282 // );
1283 // }
1284 // } else {
1285 // for(
1286 // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
1287 // mit!=iso_elastic.setOfBlocks.end();
1288 // mit++
1289 // ) {
1290 // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
1291 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1292 // ave_stress_iso.getOpPtrVector().push_back(
1293 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1294 // "DISPLACEMENT",elastic_common_data
1295 // )
1296 // );
1297 // ave_stress_iso.getOpPtrVector().push_back(
1298 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1299 // "MESH_NODE_POSITIONS",elastic_common_data
1300 // )
1301 // );
1302 // ave_stress_iso.getOpPtrVector().push_back(
1303 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1304 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,false,false,true
1305 // )
1306 // );
1307 // ave_stress_iso.getOpPtrVector().push_back(
1308 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral(
1309 // "DISPLACEMENT",elastic_block_data,elastic_common_data,stress_homo
1310 // )
1311 // );
1312 // }
1313 // for(
1314 // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
1315 // mit!=trans_elastic.setOfBlocks.end();
1316 // mit++
1317 // ) {
1318 // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
1319 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1320 // ave_stress_trans.getOpPtrVector().push_back(
1321 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1322 // "DISPLACEMENT",elastic_common_data
1323 // )
1324 // );
1325 // ave_stress_trans.getOpPtrVector().push_back(
1326 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1327 // "MESH_NODE_POSITIONS",elastic_common_data
1328 // )
1329 // );
1330 // ave_stress_trans.getOpPtrVector().push_back(
1331 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1332 // "POTENTIAL_FIELD",elastic_common_data
1333 // )
1334 // );
1335 // ave_stress_trans.getOpPtrVector().push_back(
1336 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1337 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,false,false,true
1338 // )
1339 // );
1340 // ave_stress_trans.getOpPtrVector().push_back(
1341 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral(
1342 // "DISPLACEMENT",elastic_block_data,elastic_common_data,stress_homo
1343 // )
1344 // );
1345 // }
1346 // }
1347 break;
1348 case NOHOMOBC:
1349 SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"*** Boundary conditions not set");
1350 }
1351
1352 struct MyPostProc: public PostProcVolumeOnRefinedMesh {
1353
1354 bool doPreProcess;
1355 bool doPostProcess;
1356
1357 MyPostProc(MoFEM::Interface &m_field):
1358 PostProcVolumeOnRefinedMesh(m_field),
1359 doPreProcess(true),
1360 doPostProcess(true)
1361 {}
1362
1363 void setDoPreProcess() { doPreProcess = true; }
1364 void unSetDoPreProcess() { doPreProcess = false; }
1365 void setDoPostProcess() { doPostProcess = true; }
1366 void unSetDoPostProcess() { doPostProcess = false; }
1367
1368
1369
1370 PetscErrorCode preProcess() {
1371 PetscFunctionBegin;
1372 if(doPreProcess) {
1373 ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
1374 }
1375 PetscFunctionReturn(0);
1376 }
1377 PetscErrorCode postProcess() {
1378 PetscFunctionBegin;
1379 if(doPostProcess) {
1380 ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
1381 }
1382 PetscFunctionReturn(0);
1383 }
1384 };
1385
1386 MyPostProc post_proc(m_field);
1387
1388 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
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);
1395 }
1396 for(
1397 map<int,NonlinearElasticElement::BlockData>::iterator sit = iso_elastic.setOfBlocks.begin();
1398 sit != iso_elastic.setOfBlocks.end(); sit++
1399 ) {
1400 post_proc.getOpPtrVector().push_back(
1401 new PostPorcStress(
1402 post_proc.postProcMesh,
1403 post_proc.mapGaussPts,
1404 "DISPLACEMENT",
1405 sit->second,
1406 post_proc.commonData,
1407 false
1408 )
1409 );
1410 }
1411 for(
1412 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
1413 sit != trans_elastic.setOfBlocks.end(); sit++
1414 ) {
1415 post_proc.getOpPtrVector().push_back(
1416 new PostPorcStress(
1417 post_proc.postProcMesh,
1418 post_proc.mapGaussPts,
1419 "DISPLACEMENT",
1420 sit->second,
1421 post_proc.commonData
1422 )
1423 );
1424 }
1425
1426 PetscScalar *avec;
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
1435 ); CHKERRQ(ierr);
1436 post_proc.setDoPreProcess();
1437 post_proc.unSetDoPostProcess();
1438 ierr = m_field.loop_finite_elements(
1439 "ELASTIC_MECHANICS","ELASTIC",post_proc
1440 ); CHKERRQ(ierr);
1441 post_proc.unSetDoPreProcess();
1442 post_proc.setDoPostProcess();
1443
1444 if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1445 ierr = m_field.loop_finite_elements(
1446 "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",post_proc
1447 ); CHKERRQ(ierr);
1448 {
1449 ostringstream sss;
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"
1453 ); CHKERRQ_MOAB(rval);
1454 }
1455 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
1456 if(choise_value == HOMOBCDISP) {
1457 ierr = m_field.loop_finite_elements(
1458 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCStress()
1459 ); CHKERRQ(ierr);
1460 }
1461 if(choise_value == HOMOBCTRAC) {
1462 ierr = m_field.loop_finite_elements(
1463 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCStress()
1464 ); CHKERRQ(ierr);
1465 }
1466 if(choise_value == HOMOBCPERIODIC) {
1467 ierr = m_field.loop_finite_elements(
1468 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
1469 ); CHKERRQ(ierr);
1470 }
1471 // if(choise_value == NITSCHE) {
1472 // if(stress_by_boundary_integral) {
1473 // nitsche_element_iso.addToRule = 1;
1474 // nitsche_element_trans.addToRule = 1;
1475 // ierr = m_field.loop_finite_elements(
1476 // "ELASTIC_MECHANICS","ELASTIC",nitsche_element_iso
1477 // ); CHKERRQ(ierr);
1478 // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1479 // ierr = m_field.loop_finite_elements(
1480 // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",nitsche_element_trans
1481 // ); CHKERRQ(ierr);
1482 // } else {
1483 // ave_stress_iso.addToRule = 1;
1484 // ave_stress_trans.addToRule = 1;
1485 // ierr = m_field.loop_finite_elements(
1486 // "ELASTIC_MECHANICS","ELASTIC",ave_stress_iso
1487 // ); CHKERRQ(ierr);
1488 // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1489 // ierr = m_field.loop_finite_elements(
1490 // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",ave_stress_trans
1491 // ); CHKERRQ(ierr);
1492 // }
1493 // }
1494 ierr = PetscOptionsGetReal(
1495 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
1496 ); CHKERRQ(ierr);
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];
1504 }
1505 }
1506 ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
1507
1508 PetscPrintf(
1509 PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
1510 );
1511
1512 for(int ii=0; ii<6; ii++) {
1513 PetscPrintf(
1514 PETSC_COMM_WORLD,
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)
1517 );
1518 }
1519
1520 //Saving Dmat as a bindary file to use it macro-structure
1521
1522
1523 char output_file_Dmat[255];
1524 ierr = PetscOptionsGetString(PETSC_NULL,"-my_output_file_Dmat",output_file_Dmat,255,&flg); CHKERRQ(ierr);
1525 if(flg) {
1526
1527 //Reading and writing binary files
1528 if(pcomm->rank()==0){
1529 int fd;
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);
1535 }
1536
1537 // MatrixDouble Dmat1;
1538 // Dmat1.resize(6,6); Dmat1.clear();
1539 // if(pcomm->rank()==0){
1540 // int fd;
1541 // PetscViewer view_in;
1542 // PetscViewerBinaryOpen(PETSC_COMM_SELF,output_file_Dmat,FILE_MODE_READ,&view_in);
1543 // PetscViewerBinaryGetDescriptor(view_in,&fd);
1544 // PetscBinaryRead(fd,&Dmat1(0,0),36,PETSC_DOUBLE);
1545 // PetscViewerDestroy(&view_in);
1546 // cout<< "Dmat1 After Reading= "<<Dmat1<<endl;
1547 // }
1548 //
1549 }
1550
1551 //detroy matrices
1552 for(int ii = 0;ii<6;ii++) {
1553 ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
1554 }
1555 ierr = VecDestroy(&D); CHKERRQ(ierr);
1556 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
1557 ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1558 ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
1559
1560 ierr = PetscTime(&v2);CHKERRQ(ierr);
1561 ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
1562
1563 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
1564
1565 }
1566 CATCH_ERRORS;
1567
1568
1569 MoFEM::Core::Finalize();
1570}
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:206
COL
@ COL
Definition: definitions.h:136
ROW
@ ROW
Definition: definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
NOBASE
@ NOBASE
Definition: definitions.h:72
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
H1
@ H1
continuous field
Definition: definitions.h:98
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:172
SIDESET
@ SIDESET
Definition: definitions.h:160
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:46
MOFEM_IMPOSIBLE_CASE
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
order
int order
Definition: fluid_structure_eigenproblem.cpp:68
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
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.
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::CoreInterface::loop_finite_elements
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.
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:59
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:81
MoFEM::CoreInterface::modify_problem_add_finite_element
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
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
Face_CenPos_Handle_multiIndex
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
Definition: homogenisation_periodic_atom_test.cpp:96
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:35
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:78
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:153
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:163
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:193
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1096
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:214
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:183
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:32
help
static char help[]
Definition: rve_mechanical.cpp:77
homo_bc_names
const char * homo_bc_names[]
Definition: rve_mechanical.cpp:87
BCs_RVELagrange_Disp
Definition: BCs_RVELagrange_Disp.hpp:18
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
BCs_RVELagrange_Trac_Rigid_Rot
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:22
BCs_RVELagrange_Trac_Rigid_Trans
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:18
BCs_RVELagrange_Trac
Definition: BCs_RVELagrange_Trac.hpp:22
Composite_xyzcoord
Definition: homogenisation_periodic_atom_test.cpp:72
Face_CenPos_Handle
Definition: homogenisation_periodic_atom_test.cpp:61
Hooke
Hook equation.
Definition: Hooke.hpp:33
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_field_meshset
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
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.
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:92
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
MoFEM::DeprecatedCoreInterface::seed_ref_level
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
Definition: DeprecatedCoreInterface.cpp:37
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
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...
Definition: DeprecatedCoreInterface.cpp:31
MoFEM::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition: MaterialBlocks.hpp:384
MoFEM::Mat_Elastic_TransIso::data
_data_ data
Definition: MaterialBlocks.hpp:397
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:152
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:170
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:36
MyPostProc
Definition: plot_base.cpp:52
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:570
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:40
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: PostProcOnRefMesh.hpp:707
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: PostProcOnRefMesh.hpp:808
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:968
Tri_Hand_tag
Definition: homogenisation_periodic_atom_test.cpp:71
VolumeCalculation
Calculate volume.
Definition: VolumeCalculation.hpp:31
xcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:68
ycoord_tag
Definition: homogenisation_periodic_atom_test.cpp:69
zcoord_tag
Definition: homogenisation_periodic_atom_test.cpp:70

◆ tetcircumcenter_tp()

void tetcircumcenter_tp ( double  a[3],
double  b[3],
double  c[3],
double  d[3],
double  circumcenter[3],
double *  xi,
double *  eta,
double *  zeta 
)

◆ tricircumcenter3d_tp()

void tricircumcenter3d_tp ( double  a[3],
double  b[3],
double  c[3],
double  circumcenter[3],
double *  xi,
double *  eta 
)

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 77 of file rve_mechanical.cpp.

◆ homo_bc_names

const char* homo_bc_names[]
Initial value:
= {
"disp",
"periodic",
"trac",
"nitsche",
"nohomobc"
}

Definition at line 87 of file rve_mechanical.cpp.

Generated on Tue Jul 5 2022 22:52:25 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow