93 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
96 PetscBool flg = PETSC_TRUE;
99 if(flg != PETSC_TRUE) {
100 SETERRQ(PETSC_COMM_SELF,1,
"*** ERROR -my_file (MESH FILE NEEDED)");
104 if(flg != PETSC_TRUE) {
112 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
113 if(pcomm == NULL) pcomm =
new ParallelComm(&moab,PETSC_COMM_WORLD);
116 PetscLogDouble t1,t2;
117 PetscLogDouble v1,v2;
118 ierr = PetscTime(&v1); CHKERRQ(
ierr);
119 ierr = PetscGetCPUTime(&t1); CHKERRQ(
ierr);
134 ierr = m_field.get_entities_by_ref_level(bit_level0,
BitRefLevel().set(),meshset_level0); CHKERRQ(
ierr);
151 ierr = m_field.get_cubit_msId_entities_by_dimension(103,
SIDESET,2,SurfacesFaces,
true); CHKERRQ(
ierr);
152 ierr = PetscPrintf(PETSC_COMM_WORLD,
"number of SideSet 103 = %d\n",SurfacesFaces.size()); CHKERRQ(
ierr);
157 ierr = m_field.seed_ref_level_MESHSET(BoundFacesMeshset,
BitRefLevel().set()); CHKERRQ(
ierr);
180 boost::shared_ptr<Hooke<adouble> > hooke_adouble(
new Hooke<adouble>);
181 boost::shared_ptr<Hooke<double> > hooke_double(
new Hooke<double>);
183 ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(
ierr);
184 ierr = elastic.addElement(
"ELASTIC",
"DISPLACEMENT"); CHKERRQ(
ierr);
185 ierr = elastic.setOperators(
"DISPLACEMENT",
"MESH_NODE_POSITIONS",
false,
true); CHKERRQ(
ierr);
189 lagrangian_element.addLagrangiangElement(
"LAGRANGE_ELEM",
"DISPLACEMENT",
"LAGRANGE_MUL_DISP",
"MESH_NODE_POSITIONS");
211 ierr = m_field.
loop_dofs(
"MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(
ierr);
220 ierr = m_field.build_problems(); CHKERRQ(
ierr);
227 ierr = m_field.partition_problem(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
228 ierr = m_field.partition_finite_elements(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
230 ierr = m_field.partition_ghost_dofs(
"ELASTIC_MECHANICS"); CHKERRQ(
ierr);
233 ierr = m_field.print_cubit_displacement_set(); CHKERRQ(
ierr);
234 ierr = m_field.print_cubit_force_set(); CHKERRQ(
ierr);
236 ierr = m_field.print_cubit_materials_set(); CHKERRQ(
ierr);
240 ierr = m_field.VecCreateGhost(
"ELASTIC_MECHANICS",
ROW,&
F[0]); CHKERRQ(
ierr);
241 for(
int ii = 1;ii<6;ii++) {
242 ierr = VecDuplicate(
F[0],&
F[ii]); CHKERRQ(
ierr);
246 ierr = m_field.VecCreateGhost(
"ELASTIC_MECHANICS",
COL,&
D); CHKERRQ(
ierr);
249 ierr = m_field.MatCreateMPIAIJWithArrays(
"ELASTIC_MECHANICS",&Aij); CHKERRQ(
ierr);
252 ierr = VecGhostUpdateBegin(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
253 ierr = VecGhostUpdateEnd(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
254 ierr = m_field.set_global_ghost_vector(
255 "ELASTIC_MECHANICS",
ROW,
D,INSERT_VALUES,SCATTER_REVERSE
258 for(
int ii = 0;ii<6;ii++) {
259 ierr = VecZeroEntries(
F[ii]); CHKERRQ(
ierr);
260 ierr = VecGhostUpdateBegin(
F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
261 ierr = VecGhostUpdateEnd(
F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
263 ierr = MatZeroEntries(Aij); CHKERRQ(
ierr);
266 elastic.getLoopFeLhs().snes_B = Aij;
269 lagrangian_element.setRVEBCsOperators(
"DISPLACEMENT",
"LAGRANGE_MUL_DISP",
"MESH_NODE_POSITIONS",Aij,
F);
271 cout<<
"after lagrangian_element.getLoopFeRVEBCLhs "<<endl;
273 cout<<
"after lagrangian_element.getLoopFeRVEBCRhs "<<endl;
280 ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(
ierr);
281 ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(
ierr);
283 for(
int ii = 0;ii<6;ii++) {
284 ierr = VecGhostUpdateBegin(
F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(
ierr);
285 ierr = VecGhostUpdateEnd(
F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(
ierr);
286 ierr = VecAssemblyBegin(
F[ii]); CHKERRQ(
ierr);
287 ierr = VecAssemblyEnd(
F[ii]); CHKERRQ(
ierr);
298 int volume_vec_ghost[] = { 0 };
299 ierr = VecCreateGhost(
300 PETSC_COMM_WORLD,(!m_field.
get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
302 ierr = VecZeroEntries(volume_vec); CHKERRQ(
ierr);
303 vol_elem.getOpPtrVector().push_back(
new VolumeCalculation(
"DISPLACEMENT",volume_vec));
307 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(
ierr);
308 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(
ierr);
310 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(
ierr);
312 ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(
ierr);
313 ierr = PetscPrintf(PETSC_COMM_WORLD,
"RVE Volume %3.2g\n",rve_volume); CHKERRQ(
ierr);
314 ierr = VecDestroy(&volume_vec);
327 void setDoPreProcess() { doPreProcess =
true; }
328 void unSetDoPreProcess() { doPreProcess =
false; }
329 void setDoPostProcess() { doPostProcess =
true; }
330 void unSetDoPostProcess() { doPostProcess =
false; }
339 PetscFunctionReturn(0);
346 PetscFunctionReturn(0);
352 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(
ierr);
353 ierr = post_proc.addFieldValuesPostProc(
"DISPLACEMENT"); CHKERRQ(
ierr);
354 ierr = post_proc.addFieldValuesGradientPostProc(
"DISPLACEMENT"); CHKERRQ(
ierr);
355 ierr = post_proc.addFieldValuesPostProc(
"MESH_NODE_POSITIONS"); CHKERRQ(
ierr);
356 ierr = post_proc.addFieldValuesGradientPostProc(
"MESH_NODE_POSITIONS"); CHKERRQ(
ierr);
359 map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
360 sit != elastic.setOfBlocks.end(); sit++
362 post_proc.getOpPtrVector().push_back(
364 post_proc.postProcMesh,
365 post_proc.mapGaussPts,
368 post_proc.commonData,
373 ierr = VecGhostUpdateBegin(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
374 ierr = VecGhostUpdateEnd(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
375 ierr = m_field.set_global_ghost_vector(
376 "ELASTIC_MECHANICS",
ROW,
D,INSERT_VALUES,SCATTER_REVERSE
383 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(
ierr);
384 ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(
ierr);
385 ierr = KSPSetFromOptions(solver); CHKERRQ(
ierr);
386 ierr = KSPSetUp(solver); CHKERRQ(
ierr);
395 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
396 ierr = VecCreateGhost(
397 PETSC_COMM_WORLD,(!m_field.
get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
400 lagrangian_element.setRVEBCsHomoStressOperators(
"DISPLACEMENT",
"LAGRANGE_MUL_DISP",
"MESH_NODE_POSITIONS",stress_homo);
403 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(
ierr);
404 for(
int ii = 0;ii<6;ii++) {
406 ierr = KSPSolve(solver,
F[ii],
D); CHKERRQ(
ierr);
407 ierr = VecGhostUpdateBegin(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
408 ierr = VecGhostUpdateEnd(
D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
409 ierr = m_field.set_global_ghost_vector(
410 "ELASTIC_MECHANICS",
ROW,
D,INSERT_VALUES,SCATTER_REVERSE
414 post_proc.setDoPreProcess();
415 post_proc.unSetDoPostProcess();
417 "ELASTIC_MECHANICS",
"ELASTIC",post_proc
419 post_proc.unSetDoPreProcess();
420 post_proc.setDoPostProcess();
423 sss <<
"mode_" <<
"Disp" <<
"_" << ii <<
".h5m";
424 rval = post_proc.postProcMesh.write_file(
425 sss.str().c_str(),
"MOAB",
"PARALLEL=WRITE_PART"
429 ierr = VecZeroEntries(stress_homo); CHKERRQ(
ierr);
431 "ELASTIC_MECHANICS",
"LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCStress()
434 PETSC_NULL,
"-my_rve_volume",&rve_volume,PETSC_NULL
436 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(
ierr);
437 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(
ierr);
438 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(
ierr);
439 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
440 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(
ierr);
441 for(
int jj=0; jj<6; jj++) {
442 Dmat(jj,ii) = avec[jj];
446 ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(
ierr);
448 PETSC_COMM_WORLD,
"\nHomogenised Stiffens Matrix = \n\n"
451 for(
int ii=0; ii<6; ii++) {
454 "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",
455 ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
467 myfile <<
"<<<< Homonenised stress >>>>>" << endl;
469 if(pcomm->rank()==0){
470 for(
int ii=0; ii<6; ii++){
471 myfile << boost::format(
"%.3lf") %
roundn(Dmat(ii,0)) << endl;
483 for(
int ii = 0;ii<6;ii++) {
484 ierr = VecDestroy(&
F[ii]); CHKERRQ(
ierr);
487 ierr = MatDestroy(&Aij); CHKERRQ(
ierr);
488 ierr = KSPDestroy(&solver); CHKERRQ(
ierr);
489 ierr = VecDestroy(&stress_homo); CHKERRQ(
ierr);
491 ierr = PetscTime(&v2);CHKERRQ(
ierr);
492 ierr = PetscGetCPUTime(&t2);CHKERRQ(
ierr);
494 PetscSynchronizedPrintf(PETSC_COMM_WORLD,
"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);