v0.14.0
Loading...
Searching...
No Matches
Classes | Macros | Typedefs | Functions | Variables
homogenisation_periodic_atom_test.cpp File Reference
#include <MoFEM.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#include <FEMethod_LowLevelStudent.hpp>
#include <FEMethod_UpLevelStudent.hpp>
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <Hooke.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>
#include <fstream>
#include <iostream>
#include <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <VolumeCalculation.hpp>
#include <BCs_RVELagrange_Disp.hpp>
#include <BCs_RVELagrange_Trac.hpp>
#include "BCs_RVELagrange_Periodic.hpp"
#include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
#include <boost/ptr_container/ptr_map.hpp>
#include <PostProcOnRefMesh.hpp>
#include <PostProcStresses.hpp>

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
 

Macros

#define RND_EPS   1e-6
 

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
 

Functions

double roundn (double n)
 
int main (int argc, char *argv[])
 

Variables

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

Macro Definition Documentation

◆ RND_EPS

#define RND_EPS   1e-6

Definition at line 102 of file homogenisation_periodic_atom_test.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 96 of file homogenisation_periodic_atom_test.cpp.

Function Documentation

◆ main()

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

Definition at line 137 of file homogenisation_periodic_atom_test.cpp.

137 {
138
139 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
140
141 try {
142
143 moab::Core mb_instance;
144 moab::Interface& moab = mb_instance;
145 int rank;
146 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
147
148 //Reade parameters from line command
149 PetscBool flg = PETSC_TRUE;
150 char mesh_file_name[255];
151 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
152 if(flg != PETSC_TRUE) {
153 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
154 }
155 PetscInt order;
156 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
157 if(flg != PETSC_TRUE) {
158 order = 5;
159 }
160
161 //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
162 double myapplied_strain[6];
163 int nmax=6;
164 ierr = PetscOptionsGetRealArray(PETSC_NULL,"-myapplied_strain",myapplied_strain,&nmax,&flg); CHKERRQ(ierr);
165 VectorDouble applied_strain;
166 applied_strain.resize(6);
167 cblas_dcopy(6, &myapplied_strain[0], 1, &applied_strain(0), 1);
168 // cout<<"applied_strain ="<<applied_strain<<endl;
169
170
171 //Read mesh to MOAB
172 const char *option;
173 option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
174 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
175 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
176 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
177
178 //We need that for code profiling
179 PetscLogDouble t1,t2;
180 PetscLogDouble v1,v2;
181 ierr = PetscTime(&v1); CHKERRQ(ierr);
182 ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
183
184 //Create MoFEM (Joseph) database
185 MoFEM::Core core(moab);
186 MoFEM::Interface& m_field = core;
187
188 //ref meshset ref level 0
189 ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
190
191 // stl::bitset see for more details
192 BitRefLevel bit_level0;
193 bit_level0.set(0);
194 EntityHandle meshset_level0;
195 rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
196 ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
197 ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
198
199 //=======================================================================================================
200 //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
201 //=======================================================================================================
202
203 //Populating the Multi-index container with -ve triangles
204 Range SurTrisNeg;
205 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
206 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
207 Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
208 double TriCen[3], coords_Tri[9];
209 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
210 const EntityHandle* conn_face; int num_nodes_Tri;
211
212 //get nodes attached to the face
213 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
214 //get nodal coordinates
215 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
216
217 //Find triangle centriod
218 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
219 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
220 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
221
222 //round values to 3 disimal places
223 if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0; else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0; //-ve and +ve value
224 if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0; else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
225 if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0; else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
226
227 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
228 //fill the multi-index container with centriod coordinates and triangle handles
229 Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
230 // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
231 // cout<<endl<<endl;
232 }
233
234 // double aaa;
235 // aaa=0.5011;
236 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0+0.5))/1000.0<<endl<<endl<<endl<<endl;
237 // aaa=-0.5011;
238 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0-0.5))/1000.0<<endl<<endl<<endl<<endl;
239
240
241 //Populating the Multi-index container with +ve triangles
242 Range SurTrisPos;
243 ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
244 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
245 for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
246 const EntityHandle* conn_face; int num_nodes_Tri;
247
248 //get nodes attached to the face
249 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
250 //get nodal coordinates
251 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
252
253 //Find triangle centriod
254 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
255 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
256 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
257
258 //round values to 3 disimal places
259 if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0; else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0;
260 if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0; else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
261 if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0; else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
262 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
263
264 //fill the multi-index container with centriod coordinates and triangle handles
265 Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
266 }
267
268 //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
269 double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
270 typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
271 typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
272 typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
273 Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
274 Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
275 Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
276
277 //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
278 XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
279 XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
280 YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
281 YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
282 ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
283 ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
284
285// cout<<"XcoordMin "<<XcoordMin << " XcoordMax "<<XcoordMax <<endl;
286// cout<<"YcoordMin "<<YcoordMin << " YcoordMax "<<YcoordMax <<endl;
287// cout<<"ZcoordMin "<<ZcoordMin << " ZcoordMax "<<ZcoordMax <<endl;
288 double LxRVE, LyRVE, LzRVE;
289 LxRVE=XcoordMax-XcoordMin;
290 LyRVE=YcoordMax-YcoordMin;
291 LzRVE=ZcoordMax-ZcoordMin;
292
293
294 //Creating Prisims between triangles on -ve and +ve faces
295 typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
296 Tri_Hand_iterator Tri_Neg;
297 typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
298 xyzcoord_iterator Tri_Pos;
299 Range PrismRange;
300 double XPos, YPos, ZPos;
301 //int count=0;
302
303 //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
304 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
305
306 Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
307 //cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
308
309 //corresponding +ve triangle
310 if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
311 if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
312 if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
313 Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
314 // cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
315 // cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
316
317
318 //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
319 EntityHandle PrismNodes[6];
320 vector<EntityHandle> TriNodesNeg, TriNodesPos;
321 double CoordNodeNeg[9], CoordNodePos[9];
322 rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
323 rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
324 rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
325 rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
326 for(int ii=0; ii<3; ii++){
327 PrismNodes[ii]=TriNodesNeg[ii];
328 }
329 // for(int ii=0; ii<3; ii++){
330 // cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
331 // }
332 // for(int ii=0; ii<3; ii++){
333 // cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
334 // }
335
336 //Match exact nodes to each other to avoide the problem of twisted prisms
337 double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
338 for(int ii=0; ii<3; ii++){
339 if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
340 if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
341 if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
342 for(int jj=0; jj<3; jj++){
343 //Round nodal coordinates to 3 dicimal places only for comparison
344 //round values to 3 disimal places
345 if(XNodeNeg>=0) XNodeNeg=double(int(XNodeNeg*1000.0+0.5))/1000.0; else XNodeNeg=double(int(XNodeNeg*1000.0-0.5))/1000.0;
346 if(YNodeNeg>=0) YNodeNeg=double(int(YNodeNeg*1000.0+0.5))/1000.0; else YNodeNeg=double(int(YNodeNeg*1000.0-0.5))/1000.0;
347 if(ZNodeNeg>=0) ZNodeNeg=double(int(ZNodeNeg*1000.0+0.5))/1000.0; else ZNodeNeg=double(int(ZNodeNeg*1000.0-0.5))/1000.0;
348
349 XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
350 if(XNodePos>=0) XNodePos=double(int(XNodePos*1000.0+0.5))/1000.0; else XNodePos=double(int(XNodePos*1000.0-0.5))/1000.0;
351 if(YNodePos>=0) YNodePos=double(int(YNodePos*1000.0+0.5))/1000.0; else YNodePos=double(int(YNodePos*1000.0-0.5))/1000.0;
352 if(ZNodePos>=0) ZNodePos=double(int(ZNodePos*1000.0+0.5))/1000.0; else ZNodePos=double(int(ZNodePos*1000.0-0.5))/1000.0;
353
354 if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
355 PrismNodes[3+ii]=TriNodesPos[jj];
356 break;
357 }
358 }
359 }
360 //prism nodes and their coordinates
361 double CoordNodesPrisms[18];
362 rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
363 // for(int ii=0; ii<6; ii++){
364 // cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
365 // }
366 // cout<<endl<<endl;
367 //insertion of individula prism element and its addition to range PrismRange
368 EntityHandle PeriodicPrism;
369 rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
370 PrismRange.insert(PeriodicPrism);
371
372 //Adding Prisims to Element LAGRANGE_ELEM (to loop over these prisims)
373 EntityHandle PrismRangeMeshset;
374 rval = moab.create_meshset(MESHSET_SET,PrismRangeMeshset); CHKERRQ_MOAB(rval);
375 rval = moab.add_entities(PrismRangeMeshset,PrismRange); CHKERRQ_MOAB(rval);
376 ierr = m_field.seed_ref_level_3D(PrismRangeMeshset,bit_level0); CHKERRQ(ierr);
377 ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
378
379 // //to see individual prisms
380 // Range Prism1;
381 // Prism1.insert(PeriodicPrism);
382 // EntityHandle out_meshset1;
383 // rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
384 // rval = moab.add_entities(out_meshset1,Prism1); CHKERRQ_MOAB(rval);
385 // ostringstream sss;
386 // sss << "Prism" << count << ".vtk"; count++;
387 // rval = moab.write_file(sss.str().c_str(),"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
388
389 }
390
391
392 //cout<<"PrismRange "<<PrismRange<<endl;
393 // //Saving prisms in interface.vtk
394 // EntityHandle out_meshset1;
395 // rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
396 // rval = moab.add_entities(out_meshset1,PrismRange); CHKERRQ_MOAB(rval);
397 // rval = moab.write_file("Prisms.vtk","VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
398
399
400// //=======================================================================================================
401// //=======================================================================================================
402
403 //*
404 //Define problem
405
406 //Fields
407 int field_rank=3;
408 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
409 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
410 ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr); //For lagrange multipliers to control the periodic motion
411 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //To control the rigid body motion (3 Traslations and 3 rotations)
412
413 //add entitities (by tets) to the field
414 //============================================================================================================
415 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
416 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
417
418 //Adding only -ve surfaces to the field LAGRANGE_MUL_PERIODIC (in periodic boundary conditions size of C (3M/2 x 3N))
419 //to create meshset from range SurTrisNeg
420 EntityHandle SurTrisNegMeshset;
421 rval = moab.create_meshset(MESHSET_SET,SurTrisNegMeshset); CHKERRQ_MOAB(rval);
422 rval = moab.add_entities(SurTrisNegMeshset,SurTrisNeg); CHKERRQ_MOAB(rval);
423 ierr = m_field.add_ents_to_field_by_type(SurTrisNegMeshset,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
424 //============================================================================================================
425
426 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
427 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
428 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
429 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
430
431 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
432 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
433 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
434 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
435
436 //ierr = m_field.set_field_order(0,MBPRISM,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
437 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
438 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
439 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
440
441 //Define FE
442 //define eleatic element
443 boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>);
444 boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>);
445 NonlinearElasticElement elastic(m_field,2);
446 ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(ierr);
447 ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
448 ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
449
450 BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
451 BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
452
453 lagrangian_element_periodic.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",PrismRange);
454 lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
455
456 //define problems
457 ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
458
459 //set finite elements for problem
460 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
461 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
462 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
463
464
465 //set refinement level for problem
466 ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",bit_level0); CHKERRQ(ierr);
467
468 /***/
469 //Declare problem
470 /****/
471 //build database
472
473 //build field
474 ierr = m_field.build_fields(); CHKERRQ(ierr);
475 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
476 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
477
478
479 //build finite elemnts
480 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
481
482 //build adjacencies
483 ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
484
485 //build problem
486 ierr = m_field.build_problems(); CHKERRQ(ierr);
487
488 /****/
489 //mesh partitioning
490
491 //partition
492 ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
493 ierr = m_field.partition_finite_elements("ELASTIC_MECHANICS"); CHKERRQ(ierr);
494 //what are ghost nodes, see Petsc Manual
495 ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
496
497 //print bcs
498 ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
499 ierr = m_field.print_cubit_force_set(); CHKERRQ(ierr);
500 //print block sets with materials
501 ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
502
503 //create matrices (here F, D and Aij are matrices for the full problem)
504 Vec D;
505 vector<Vec> F(6);
506 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
507 for(int ii = 1;ii<6;ii++) {
508 ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
509 }
510
511 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
512 Mat Aij;
513 ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
514
515 for(int ii = 0;ii<6;ii++) {
516 ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
517 ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
518 ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
519 }
520
521 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
522
523
524 //Assemble F and Aij
525 elastic.getLoopFeLhs().snes_B = Aij;
526 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs()); CHKERRQ(ierr);
527
528 lagrangian_element_periodic.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,F);
529 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
530 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
531
532 BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
533 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element_periodic.setOfRVEBC);
534 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
535
536 for(int ii = 0;ii<6;ii++) {
537 ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
538 ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
539 ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
540 ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
541 }
542
543 ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
544 ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
545
546 // Volume calculation
547 //==============================================================================================================================
548 VolumeElementForcesAndSourcesCore vol_elem(m_field);
549 Vec volume_vec;
550 int volume_vec_ghost[] = { 0 };
551 ierr = VecCreateGhost(
552 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
553 ); CHKERRQ(ierr);
554 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
555 vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
556
557 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC", vol_elem); CHKERRQ(ierr);
558
559 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
560 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
561 double rve_volume;
562 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
563
564 ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
565 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
566 ierr = VecDestroy(&volume_vec);
567
568
569 //==============================================================================================================================
570 //Post processing
571 //==============================================================================================================================
572 struct MyPostProc: public PostProcVolumeOnRefinedMesh {
573 bool doPreProcess;
574 bool doPostProcess;
575 MyPostProc(MoFEM::Interface &m_field):
576 PostProcVolumeOnRefinedMesh(m_field),
577 doPreProcess(true),
578 doPostProcess(true)
579 {}
580 void setDoPreProcess() { doPreProcess = true; }
581 void unSetDoPreProcess() { doPreProcess = false; }
582 void setDoPostProcess() { doPostProcess = true; }
583 void unSetDoPostProcess() { doPostProcess = false; }
584
585
586
587 PetscErrorCode preProcess() {
588 PetscFunctionBegin;
589 if(doPreProcess) {
590 ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
591 }
592 PetscFunctionReturn(0);
593 }
594 PetscErrorCode postProcess() {
595 PetscFunctionBegin;
596 if(doPostProcess) {
597 ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
598 }
599 PetscFunctionReturn(0);
600 }
601 };
602
603 MyPostProc post_proc(m_field);
604
605 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
606 ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
607 ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
608 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
609 ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
610
611 for(
612 map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
613 sit != elastic.setOfBlocks.end(); sit++
614 ) {
615 post_proc.getOpPtrVector().push_back(
616 new PostPorcStress(
617 post_proc.postProcMesh,
618 post_proc.mapGaussPts,
619 "DISPLACEMENT",
620 sit->second,
621 post_proc.commonData,
622 false
623 )
624 );
625 }
626 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
627 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
628 ierr = m_field.set_global_ghost_vector(
629 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
630 ); CHKERRQ(ierr);
631
632 //==============================================================================================================================
633
634 //Solver
635 KSP solver;
636 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
637 ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
638 ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
639 ierr = KSPSetUp(solver); CHKERRQ(ierr);
640
641 MatrixDouble Dmat;
642 Dmat.resize(6,6);
643 Dmat.clear();
644
645 //calculate honmogenised stress
646 //create a vector for 6 components of homogenized stress
647 Vec stress_homo;
648 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
649 ierr = VecCreateGhost(
650 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
651 ); CHKERRQ(ierr);
652
653 lagrangian_element_periodic.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo);
654
655 PetscScalar *avec;
656 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
657 for(int ii = 0;ii<6;ii++) {
658 ierr = VecZeroEntries(D); CHKERRQ(ierr);
659 ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
660 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
661 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
662 ierr = m_field.set_global_ghost_vector(
663 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
664 ); CHKERRQ(ierr);
665
666 //post-processing the resutls
667 post_proc.setDoPreProcess();
668 post_proc.unSetDoPostProcess();
669 ierr = m_field.loop_finite_elements(
670 "ELASTIC_MECHANICS","ELASTIC",post_proc
671 ); CHKERRQ(ierr);
672 post_proc.unSetDoPreProcess();
673 post_proc.setDoPostProcess();
674 {
675 ostringstream sss;
676 sss << "mode_" << "Disp" << "_" << ii << ".h5m";
677 rval = post_proc.postProcMesh.write_file(
678 sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
679 ); CHKERRQ_MOAB(rval);
680 }
681
682 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
683 ierr = m_field.loop_finite_elements(
684 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
685 ); CHKERRQ(ierr);
686 ierr = PetscOptionsGetReal(
687 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
688 ); CHKERRQ(ierr);
689 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
690 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
691 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
692 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
693 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
694 for(int jj=0; jj<6; jj++) {
695 Dmat(jj,ii) = avec[jj];
696 }
697 }
698
699 ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
700 PetscPrintf(
701 PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
702 );
703
704 for(int ii=0; ii<6; ii++) {
705 PetscPrintf(
706 PETSC_COMM_WORLD,
707 "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",
708 ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
709 );
710 }
711 //==============================================================================================================================
712
713 //Open mesh_file_name.txt for writing
714 ofstream myfile;
715 myfile.open ((string(mesh_file_name)+".txt").c_str());
716
717 //Output displacements
718 myfile << "<<<< Homonenised stress >>>>>" << endl;
719
720 if(pcomm->rank()==0){
721 for(int ii=0; ii<6; ii++){
722 myfile << boost::format("%.3lf") % roundn(Dmat(ii,0)) << endl;
723 avec++;
724 }
725 }
726 //Close mesh_file_name.txt
727 myfile.close();
728
729
730
731 //detroy matrices
732 for(int ii = 0;ii<6;ii++) {
733 ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
734 }
735 ierr = VecDestroy(&D); CHKERRQ(ierr);
736 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
737 ierr = KSPDestroy(&solver); CHKERRQ(ierr);
738 ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
739
740 ierr = PetscTime(&v2);CHKERRQ(ierr);
741 ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
742
743 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
744
745 }
746 CATCH_ERRORS;
747
748 MoFEM::Core::Finalize();
749
750}
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
SIDESET
@ SIDESET
Definition definitions.h:160
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition definitions.h:467
order
constexpr int order
Definition dg_projection.cpp:18
F
@ F
Definition free_surface.cpp:394
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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::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.
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
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
roundn
double roundn(double n)
Definition homogenisation_periodic_atom_test.cpp:105
help
static char help[]
Definition homogenisation_periodic_atom_test.cpp:54
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::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::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition DeprecatedPetsc.hpp:192
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
BCs_RVELagrange_Periodic
Definition BCs_RVELagrange_Periodic.hpp:19
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::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_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::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
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition HookeElement.hpp:27
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

◆ roundn()

double roundn ( double n)

Definition at line 105 of file homogenisation_periodic_atom_test.cpp.

106{
107 //break n into fractional part (fract) and integral part (intp)
108 double fract, intp;
109 fract = modf(n,&intp);
110
111// //round up
112// if (fract>=.5)
113// {
114// n*=10;
115// ceil(n);
116// n/=10;
117// }
118// //round down
119// if (fract<=.5)
120// {
121// n*=10;
122// floor(n);
123// n/=10;
124// }
125 // case where n approximates zero, set n to "positive" zero
126 if (abs(intp)==0)
127 {
128 if(abs(fract)<=RND_EPS)
129 {
130 n=0.000;
131 }
132 }
133 return n;
134}
RND_EPS
#define RND_EPS
Definition homogenisation_periodic_atom_test.cpp:102
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38

Variable Documentation

◆ help

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

Definition at line 54 of file homogenisation_periodic_atom_test.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow