v0.14.0
Loading...
Searching...
No Matches
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 <NitscheMethod.hpp>
#include <moab/AdaptiveKDTree.hpp>
#include <NitschePeriodicMethod.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:193
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 nme:nme847.
Definition definitions.h:60
NOBASE
@ NOBASE
Definition definitions.h:59
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:84
H1
@ H1
continuous field
Definition definitions.h:85
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:33
order
constexpr int order
Definition dg_projection.cpp:18
F
@ F
Definition free_surface.cpp:394
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
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::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::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_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::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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:49
_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:71
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
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
_IT_NUMEREDDOF_ROW_BY_ENT_FOR_LOOP_
#define _IT_NUMEREDDOF_ROW_BY_ENT_FOR_LOOP_(PROBLEMPTR, ENT, IT)
use with loops to iterate row DOFs
Definition ProblemsMultiIndices.hpp:281
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
D
double D
Definition initial_diffusion.cpp:44
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:152
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition DeprecatedPetsc.hpp:203
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition EntitiesFieldData.hpp:1589
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:21
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:82
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:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
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:24
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:18
MoFEM::ForcesAndSourcesCore::preProcess
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition ForcesAndSourcesCore.cpp:2011
MoFEM::ForcesAndSourcesCore::postProcess
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
Definition ForcesAndSourcesCore.cpp:2034
MoFEM::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition MaterialBlocks.hpp:371
MoFEM::Mat_Elastic_TransIso::data
_data_ data
Definition MaterialBlocks.hpp:384
MoFEM::Mat_Elastic
Elastic material data structure.
Definition MaterialBlocks.hpp:139
MoFEM::Mat_Elastic::data
_data_ data
Definition MaterialBlocks.hpp:157
MoFEM::Problem
keeps basic data about problem
Definition ProblemsMultiIndices.hpp:54
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
MyPostProc
Definition plot_base.cpp:34
MyPostProc::postProcess
MoFEMErrorCode postProcess()
Definition plot_base.cpp:651
MyPostProc::preProcess
MoFEMErrorCode preProcess()
Definition plot_base.cpp:642
NitscheMethod::BlockData
Block data for Nitsche method.
Definition NitscheMethod.hpp:82
NitscheMethod::BlockData::gamma
double gamma
Penalty term, see nitsche_method_hal.
Definition NitscheMethod.hpp:83
NitscheMethod::BlockData::faceElemName
string faceElemName
name of element face
Definition NitscheMethod.hpp:85
NitscheMethod::BlockData::phi
double phi
Nitsche method parameter, see nitsche_method_hal.
Definition NitscheMethod.hpp:84
NitscheMethod::BlockData::fAces
Range fAces
faces on which constrain is applied
Definition NitscheMethod.hpp:86
NitscheMethod::CommonData
Common data shared between finite element operators.
Definition NitscheMethod.hpp:91
NitscheMethod::OpLhsNormal
Calculate Nitsche method terms on left hand side.
Definition NitscheMethod.hpp:513
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition NonLinearElasticElement.hpp:84
NonlinearElasticElement::CommonData
common data used by volume elements
Definition NonLinearElasticElement.hpp:105
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition NonLinearElasticElement.hpp:557
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition HookeElement.hpp:27
PeriodicNitscheConstrains::OpRhsPeriodicNormal::calcualteDMatrix
static PetscErrorCode calcualteDMatrix(VectorAdaptor coords, MatrixDouble &mat_d)
Definition NitschePeriodicMethod.hpp:612
PostProcVolumeOnRefinedMesh
Post processing.
Definition PostProcOnRefMesh.hpp:957
Tri_Hand_tag
Definition homogenisation_periodic_atom_test.cpp:71
VolumeCalculation
Calculate volume.
Definition VolumeCalculation.hpp:18
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.

87 {
88 "disp",
89 "periodic",
90 "trac",
91 "nitsche",
92 "nohomobc"
93};
Generated by Doxygen 1.12.0 and hosted at University of Glasgow