v0.15.0
Loading...
Searching...
No Matches
Macros | Functions | Variables
homogenisation_trac_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_Trac_Rigid_Rot.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.

Macros

#define RND_EPS   1e-6
 

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 52 of file homogenisation_trac_atom_test.cpp.

Function Documentation

◆ main()

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

Definition at line 86 of file homogenisation_trac_atom_test.cpp.

86 {
87
88 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
89
90 try {
91
92 moab::Core mb_instance;
93 moab::Interface& moab = mb_instance;
94 int rank;
95 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
96
97 //Reade parameters from line command
98 PetscBool flg = PETSC_TRUE;
99 char mesh_file_name[255];
100 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
101 if(flg != PETSC_TRUE) {
102 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
103 }
104 PetscInt order;
105 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
106 if(flg != PETSC_TRUE) {
107 order = 5;
108 }
109
110 //Read mesh to MOAB
111 const char *option;
112 option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
113 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
114 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
115 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
116
117 //We need that for code profiling
118 PetscLogDouble t1,t2;
119 PetscLogDouble v1,v2;
120 ierr = PetscTime(&v1); CHKERRQ(ierr);
121 ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
122
123 //Create MoFEM (Joseph) database
124 MoFEM::Core core(moab);
125 MoFEM::Interface& m_field = core;
126
127 //ref meshset ref level 0
128 ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
129
130 // stl::bitset see for more details
131 BitRefLevel bit_level0;
132 bit_level0.set(0);
133 EntityHandle meshset_level0;
134 rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
135 ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
136 ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
137
138 /***/
139 //Define problem
140
141 //Fields
142 int field_rank=3;
143 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
144 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
145 ierr = m_field.add_field("LAGRANGE_MUL_TRAC",NOFIELD,NOBASE,6); CHKERRQ(ierr);
146 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //Control 3 rigid body translations in x, y and z axis
147 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_ROT",NOFIELD,NOBASE,3); CHKERRQ(ierr); //Controla 3 rigid body rotations about x, y and z axis
148
149
150 //add entitities (by tets) to the field
151 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
152 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
153
154 //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
155 //int order = 5;
156 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
157 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
158 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
159 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
160
161 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
162 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
163 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
164 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
165
166
167 //FE
168 //Define FE
169 //define eleatic element
170 boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>);
171 boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>);
172 NonlinearElasticElement elastic(m_field,2);
173 ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(ierr);
174 ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
175 ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
176
177
178 BCs_RVELagrange_Trac lagrangian_element(m_field);
179 lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS");
180 lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
181 lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM_ROT","DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT","MESH_NODE_POSITIONS");
182
183 //define problems
184 ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
185
186 //set finite elements for problem
187 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
188 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
189 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
190 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT"); CHKERRQ(ierr);
191
192
193 //set refinement level for problem
194 ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",bit_level0); CHKERRQ(ierr);
195
196 /***/
197 //Declare problem
198 /****/
199 //build database
200
201 //build field
202 ierr = m_field.build_fields(); CHKERRQ(ierr);
203 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
204 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
205
206 //build finite elemnts
207 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
208
209 //build adjacencies
210 ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
211
212 //build problem
213 ierr = m_field.build_problems(); CHKERRQ(ierr);
214
215 /****/
216 //mesh partitioning
217
218 //partition
219 ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
220 ierr = m_field.partition_finite_elements("ELASTIC_MECHANICS"); CHKERRQ(ierr);
221 //what are ghost nodes, see Petsc Manual
222 ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
223
224 //print bcs
225 ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
226 ierr = m_field.print_cubit_force_set(); CHKERRQ(ierr);
227 //print block sets with materials
228 ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
229
230 //create matrices (here F, D and Aij are matrices for the full problem)
231 vector<Vec> F(6);
232 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
233 for(int ii = 1;ii<6;ii++) {
234 ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
235 }
236
237 Vec D;
238 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
239
240 Mat Aij;
241 ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
242
243 ierr = VecZeroEntries(D); CHKERRQ(ierr);
244 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
245 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
246 ierr = m_field.set_global_ghost_vector(
247 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
248 ); CHKERRQ(ierr);
249
250 for(int ii = 0;ii<6;ii++) {
251 ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
252 ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
253 ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
254 }
255 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
256
257 //Assemble F and Aij
258 elastic.getLoopFeLhs().snes_B = Aij;
259 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs()); CHKERRQ(ierr);
260 lagrangian_element.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",Aij,F);
261 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
262 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
263
264 BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
265 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element.setOfRVEBC);
266 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
267
268 BCs_RVELagrange_Trac_Rigid_Rot lagrangian_element_rigid_body_rot(m_field);
269 lagrangian_element_rigid_body_rot.setRVEBCsRigidBodyRotOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT",Aij, lagrangian_element.setOfRVEBC);
270 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT",lagrangian_element_rigid_body_rot.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
271
272// //Matrix View
273// MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
274// std::string wait;
275// std::cin >> wait;
276
277 ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
278 ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
279
280 for(int ii = 0;ii<6;ii++) {
281 ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
282 ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
283 ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
284 ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
285 }
286
287// //ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
288// //ierr = MatView(Aij,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
289//
290////Matrix View
291//MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
292//std::string wait;
293//std::cin >> wait;
294
295
296
297 // Volume calculation
298 //==============================================================================================================================
299 VolumeElementForcesAndSourcesCore vol_elem(m_field);
300 Vec volume_vec;
301 int volume_vec_ghost[] = { 0 };
302 ierr = VecCreateGhost(
303 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
304 ); CHKERRQ(ierr);
305 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
306 vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
307
308 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC", vol_elem); CHKERRQ(ierr);
309
310 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
311 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
312 double rve_volume;
313 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
314
315 ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
316 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
317 ierr = VecDestroy(&volume_vec);
318
319 //==============================================================================================================================
320 //Post processing
321 //==============================================================================================================================
322 struct MyPostProc: public PostProcVolumeOnRefinedMesh {
323 bool doPreProcess;
324 bool doPostProcess;
325 MyPostProc(MoFEM::Interface &m_field):
326 PostProcVolumeOnRefinedMesh(m_field),
327 doPreProcess(true),
328 doPostProcess(true)
329 {}
330 void setDoPreProcess() { doPreProcess = true; }
331 void unSetDoPreProcess() { doPreProcess = false; }
332 void setDoPostProcess() { doPostProcess = true; }
333 void unSetDoPostProcess() { doPostProcess = false; }
334
335
336
337 PetscErrorCode preProcess() {
338 PetscFunctionBegin;
339 if(doPreProcess) {
340 ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
341 }
342 PetscFunctionReturn(0);
343 }
344 PetscErrorCode postProcess() {
345 PetscFunctionBegin;
346 if(doPostProcess) {
347 ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
348 }
349 PetscFunctionReturn(0);
350 }
351 };
352
353 MyPostProc post_proc(m_field);
354
355 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
356 ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
357 ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
358 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
359 ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
360
361 for(
362 map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
363 sit != elastic.setOfBlocks.end(); sit++
364 ) {
365 post_proc.getOpPtrVector().push_back(
366 new PostPorcStress(
367 post_proc.postProcMesh,
368 post_proc.mapGaussPts,
369 "DISPLACEMENT",
370 sit->second,
371 post_proc.commonData,
372 false
373 )
374 );
375 }
376 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
377 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
378 ierr = m_field.set_global_ghost_vector(
379 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
380 ); CHKERRQ(ierr);
381
382 //==============================================================================================================================
383
384 //Solver
385 KSP solver;
386 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
387 ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
388 ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
389 ierr = KSPSetUp(solver); CHKERRQ(ierr);
390
391 MatrixDouble Dmat;
392 Dmat.resize(6,6);
393 Dmat.clear();
394
395 //calculate honmogenised stress
396 //create a vector for 6 components of homogenized stress
397 Vec stress_homo;
398 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
399 ierr = VecCreateGhost(
400 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
401 ); CHKERRQ(ierr);
402
403 lagrangian_element.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",stress_homo);
404
405 PetscScalar *avec;
406 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
407 for(int ii = 0;ii<6;ii++) {
408 ierr = VecZeroEntries(D); CHKERRQ(ierr);
409 ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
410 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
411 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
412 ierr = m_field.set_global_ghost_vector(
413 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
414 ); CHKERRQ(ierr);
415
416 //post-processing the resutls
417 post_proc.setDoPreProcess();
418 post_proc.unSetDoPostProcess();
419 ierr = m_field.loop_finite_elements(
420 "ELASTIC_MECHANICS","ELASTIC",post_proc
421 ); CHKERRQ(ierr);
422 post_proc.unSetDoPreProcess();
423 post_proc.setDoPostProcess();
424 {
425 ostringstream sss;
426 sss << "mode_" << "Disp" << "_" << ii << ".h5m";
427 rval = post_proc.postProcMesh.write_file(
428 sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
429 ); CHKERRQ_MOAB(rval);
430 }
431
432 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
433 ierr = m_field.loop_finite_elements(
434 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCStress()
435 ); CHKERRQ(ierr);
436 ierr = PetscOptionsGetReal(
437 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
438 ); CHKERRQ(ierr);
439 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
440 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
441 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
442 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
443 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
444 for(int jj=0; jj<6; jj++) {
445 Dmat(jj,ii) = avec[jj];
446 }
447 }
448
449 ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
450 PetscPrintf(
451 PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
452 );
453
454 for(int ii=0; ii<6; ii++) {
455 PetscPrintf(
456 PETSC_COMM_WORLD,
457 "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",
458 ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
459 );
460 }
461 //==============================================================================================================================
462
463 //Open mesh_file_name.txt for writing
464 ofstream myfile;
465 myfile.open ((string(mesh_file_name)+".txt").c_str());
466
467 //Output displacements
468 myfile << "<<<< Homonenised stress >>>>>" << endl;
469
470 if(pcomm->rank()==0){
471 for(int ii=0; ii<6; ii++){
472 myfile << boost::format("%.3lf") % roundn(Dmat(ii,0)) << endl;
473 avec++;
474 }
475 }
476 //Close mesh_file_name.txt
477 myfile.close();
478
479
480
481 //detroy matrices
482 for(int ii = 0;ii<6;ii++) {
483 ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
484 }
485 ierr = VecDestroy(&D); CHKERRQ(ierr);
486 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
487 ierr = KSPDestroy(&solver); CHKERRQ(ierr);
488 ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
489
490 ierr = PetscTime(&v2);CHKERRQ(ierr);
491 ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
492
493 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
494
495 }
496 CATCH_ERRORS;
497
498 MoFEM::Core::Finalize();
499
500}
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 .
Definition definitions.h:60
NOBASE
@ NOBASE
Definition definitions.h:59
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
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition definitions.h:470
order
constexpr int order
Order displacement.
Definition dg_projection.cpp:38
F
@ F
Definition free_surface.cpp:624
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
roundn
double roundn(double n)
Definition homogenisation_trac_atom_test.cpp:55
help
static char help[]
Definition homogenisation_trac_atom_test.cpp:51
D
double D
Definition initial_diffusion.cpp:44
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition MatrixFunction.hpp:64
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::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
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
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:118
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::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
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
VolumeCalculation
Calculate volume.
Definition VolumeCalculation.hpp:18

◆ roundn()

double roundn ( double  n)

Definition at line 55 of file homogenisation_trac_atom_test.cpp.

56{
57 //break n into fractional part (fract) and integral part (intp)
58 double fract, intp;
59 fract = modf(n,&intp);
60
61// //round up
62// if (fract>=.5)
63// {
64// n*=10;
65// ceil(n);
66// n/=10;
67// }
68// //round down
69// if (fract<=.5)
70// {
71// n*=10;
72// floor(n);
73// n/=10;
74// }
75 // case where n approximates zero, set n to "positive" zero
76 if (abs(intp)==0)
77 {
78 if(abs(fract)<=RND_EPS)
79 {
80 n=0.000;
81 }
82 }
83 return n;
84}
RND_EPS
#define RND_EPS
Definition homogenisation_trac_atom_test.cpp:52
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 51 of file homogenisation_trac_atom_test.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow