v0.14.0
Loading...
Searching...
No Matches
Classes | Typedefs | Functions | Variables
rve_mech_plas_interface_periodic.cpp File Reference
#include <MoFEM.hpp>
#include <boost/program_options.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#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 <CohesiveInterfaceElement.hpp>
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <SmallTransverselyIsotropic.hpp>
#include <SmallStrainPlasticity.hpp>
#include <SmallStrainPlasticityMaterialModels.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <PostProcOnRefMesh.hpp>
#include <string>

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
 

Functions

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

Variables

static char help []
 

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 115 of file rve_mech_plas_interface_periodic.cpp.

Function Documentation

◆ main()

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

Definition at line 122 of file rve_mech_plas_interface_periodic.cpp.

122 {
123 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
124
125 try {
126
127 moab::Core mb_instance;
128 moab::Interface& moab = mb_instance;
129 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
130 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
131
132 PetscBool flg = PETSC_TRUE;
133 char mesh_file_name[255];
134 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
135 if(flg != PETSC_TRUE) {
136 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
137 }
138
139 PetscInt order;
140 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
141 if(flg != PETSC_TRUE) {
142 order = 2;
143 }
144 PetscInt bubble_order;
145 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_bubble_order",&bubble_order,&flg); CHKERRQ(ierr);
146 if(flg != PETSC_TRUE) {
147 bubble_order = order;
148 }
149
150 // use this if your mesh is partitioned and you run code on parts,
151 // you can solve very big problems
152 PetscBool is_partitioned = PETSC_FALSE;
153 ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
154
155 //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
156 double mygiven_strain[6];
157 int nmax=6;
158 ierr = PetscOptionsGetRealArray(PETSC_NULL,"-my_given_strain",mygiven_strain,&nmax,&flg); CHKERRQ(ierr);
159 VectorDouble given_strain;
160 given_strain.resize(6);
161 cblas_dcopy(6, &mygiven_strain[0], 1, &given_strain(0), 1);
162 cout<<"given_strain ="<<given_strain<<endl;
163
164 //Read mesh to MOAB
165
166 if(is_partitioned == PETSC_TRUE) {
167 const char *option;
168 option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
169 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
170 } else {
171 const char *option;
172 option = "";
173 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
174 }
175
176
177
178 double young_modulus = 1;
179 SmallStrainParaboloidalPlasticity cp;
180 {
181 cp.tAgs.resize(3);
182 cp.tAgs[0] = 3;
183 cp.tAgs[1] = 4;
184 cp.tAgs[2] = 5;
185 cp.tOl = 1e-12;
186
187 double poisson_ratio = 0.25;
188 cp.sIgma_yt = 1;
189 cp.sIgma_yc = 1;
190
191 cp.Ht = 0.1;
192 cp.Hc = 0.1;
193
194 cp.nup = 0.3;
195
196
197 {
198 ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
199 ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
200 cp.mu = MU(young_modulus,poisson_ratio);
201 cp.lambda = LAMBDA(young_modulus,poisson_ratio);
202 ierr = PetscOptionsGetReal(0,"-my_sigma_yt",&cp.sIgma_yt,0); CHKERRQ(ierr);
203 ierr = PetscOptionsGetReal(0,"-my_sigma_yc",&cp.sIgma_yc,0); CHKERRQ(ierr);
204
205 ierr = PetscOptionsGetReal(0,"-my_Ht",&cp.Ht,0); CHKERRQ(ierr);
206 ierr = PetscOptionsGetReal(0,"-my_Hc",&cp.Hc,0); CHKERRQ(ierr);
207
208 ierr = PetscOptionsGetReal(0,"-my_nt",&cp.nt,0); CHKERRQ(ierr);
209 ierr = PetscOptionsGetReal(0,"-my_nc",&cp.nc,0); CHKERRQ(ierr);
210
211
212 ierr = PetscOptionsGetReal(0,"-my_nup",&cp.nup,0); CHKERRQ(ierr);
213 }
214
215 PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
216 PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
217
218 PetscPrintf(PETSC_COMM_WORLD,"sIgma_yt = %4.2e \n",cp.sIgma_yt);
219 PetscPrintf(PETSC_COMM_WORLD,"sIgma_yc = %4.2e \n",cp.sIgma_yc);
220
221 PetscPrintf(PETSC_COMM_WORLD,"nup = %4.2e \n",cp.nup);
222
223 cp.sTrain.resize(6,false);
224 cp.sTrain.clear();
225 cp.plasticStrain.resize(6,false);
226 cp.plasticStrain.clear();
227 cp.internalVariables.resize(2,false);
228 cp.internalVariables.clear();
229 cp.createMatAVecR();
230 cp.snesCreate();
231 // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
232 }
233
234
235 // double young_modulus = 1;
236 // SmallStrainJ2Plasticity cp;
237 // {
238 // cp.tAgs.resize(3);
239 // cp.tAgs[0] = 3;
240 // cp.tAgs[1] = 4;
241 // cp.tAgs[2] = 5;
242 // cp.tOl = 1e-12;
243 //
244 // double poisson_ratio = 0.25;
245 // cp.sIgma_y = 1;
246 // cp.H = 0.1;
247 // cp.K = 0;
248 // cp.phi = 1;
249 // {
250 // ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
251 // ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
252 // cp.mu = MU(young_modulus,poisson_ratio);
253 // cp.lambda = LAMBDA(young_modulus,poisson_ratio);
254 // ierr = PetscOptionsGetReal(0,"-my_sigma_y",&cp.sIgma_y,0); CHKERRQ(ierr);
255 // ierr = PetscOptionsGetReal(0,"-my_H",&cp.H,0); CHKERRQ(ierr);
256 // ierr = PetscOptionsGetReal(0,"-my_K",&cp.K,0); CHKERRQ(ierr);
257 // ierr = PetscOptionsGetReal(0,"-my_phi",&cp.phi,0); CHKERRQ(ierr);
258 // }
259 //
260 // PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
261 // PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
262 // PetscPrintf(PETSC_COMM_WORLD,"sIgma_y = %4.2e \n",cp.sIgma_y);
263 // PetscPrintf(PETSC_COMM_WORLD,"H = %4.2e \n",cp.H);
264 // PetscPrintf(PETSC_COMM_WORLD,"K = %4.2e \n",cp.K);
265 // PetscPrintf(PETSC_COMM_WORLD,"phi = %4.2e \n",cp.phi);
266 //
267 //
268 // cp.sTrain.resize(6,false);
269 // cp.sTrain.clear();
270 // cp.plasticStrain.resize(6,false);
271 // cp.plasticStrain.clear();
272 // cp.internalVariables.resize(7,false);
273 // cp.internalVariables.clear();
274 // cp.createMatAVecR();
275 // cp.snesCreate();
276 // // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
277 // }
278
279
280 //Create MoFEM (Joseph) database
281 MoFEM::Core core(moab);
282 MoFEM::Interface& m_field = core;
283
284 vector<BitRefLevel> bit_levels;
285 Tag th_meshset_info;
286 int def_meshset_info[2] = {0,0};
287 rval = moab.tag_get_handle("MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info);
288 int meshset_data[2];
289 EntityHandle root = moab.get_root_set();
290 rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
291 if(meshset_data[0]==0) {
292 meshset_data[0] = 1;
293 rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
294 }
295 bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
296
297
298// EntityHandle out_meshset;
299// rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
300// // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
301// ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
302// Range LatestRefinedTets;
303// rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
304// Range LatestRefinedPrisms;
305// rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
306
307
308 BitRefLevel problem_bit_level = bit_levels.back();
309
310 // const clock_t begin_time = clock();
311 ierr = m_field.build_fields(); CHKERRQ(ierr);
312 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
313
314
315 Range interface_prims_on_problem_bit_level;
316 ierr = m_field.get_entities_by_type_and_ref_level(
317 problem_bit_level,BitRefLevel().set(),MBPRISM,interface_prims_on_problem_bit_level
318 ); CHKERRQ(ierr);
319 Range tets_on_problem_bit_level;
320 ierr = m_field.get_entities_by_type_and_ref_level(
321 problem_bit_level,BitRefLevel().set(),MBTET,tets_on_problem_bit_level
322 ); CHKERRQ(ierr);
323
324 cout<<"interface_prims_on_problem_bit_level "<<interface_prims_on_problem_bit_level.size()<<endl;
325 cout<<"tets_on_problem_bit_level "<<tets_on_problem_bit_level.size()<<endl;
326
327 //=======================================================================================================
328 //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
329 //=======================================================================================================
330 Range preriodic_prisms;
331 //Populating the Multi-index container with -ve triangles
332
333
334 Range SurTrisNeg; //Old + new traiangles on side set 101
335 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
336 ierr = PetscPrintf(PETSC_COMM_WORLD,"Old+New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
337
338 Range AllTris; //All new traiangles in the new-bitlevel
339 ierr = m_field.get_entities_by_type_and_ref_level(problem_bit_level,BitRefLevel().set(),MBTRI,AllTris); CHKERRQ(ierr);
340 ierr = PetscPrintf(PETSC_COMM_WORLD,"All triangles in the new bit-level = %d\n",AllTris.size()); CHKERRQ(ierr);
341
342 //Only triangles in the new bit-level
343 SurTrisNeg = intersect(AllTris,SurTrisNeg);
344 ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
345
346// cout<<"SurTrisNeg "<<SurTrisNeg<<endl;
347 Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
348 double TriCen[3], coords_Tri[9];
349
350
351
352 double roundfact=1e3;
353
354 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
355 const EntityHandle* conn_face; int num_nodes_Tri;
356
357 //get nodes attached to the face
358 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
359 //get nodal coordinates
360 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
361
362 //Find triangle centriod
363 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
364 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
365 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
366
367// cout<<"\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
368 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
369 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
370 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
371
372// cout<<"TriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
373 //fill the multi-index container with centriod coordinates and triangle handles
374 Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
375 // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
376 // cout<<endl<<endl;
377 }
378
379 // double aaa;
380 // aaa=0.5011;
381 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*roundfact+0.5))/roundfact<<endl<<endl<<endl<<endl;
382 // aaa=-0.5011;
383 // cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*roundfact-0.5))/roundfact<<endl<<endl<<endl<<endl;
384
385
386 //Populating the Multi-index container with +ve triangles
387 Range SurTrisPos;
388 ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
389 ierr = PetscPrintf(PETSC_COMM_WORLD,"Old+New triangles on SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
390
391 //Only triangles in the new bit-level
392 SurTrisPos = intersect(AllTris,SurTrisPos);
393 ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
394
395
396
397
398 for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
399 const EntityHandle* conn_face; int num_nodes_Tri;
400
401 //get nodes attached to the face
402 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
403 //get nodal coordinates
404 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
405
406 //Find triangle centriod
407 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
408 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
409 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
410
411 //round values to 3 disimal places
412 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
413 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
414 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
415
416// cout<<"TriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
417
418 //fill the multi-index container with centriod coordinates and triangle handles
419 Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
420 }
421
422 //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
423 double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
424 typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
425 typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
426 typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
427 Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
428 Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
429 Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
430
431 //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
432 XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
433 XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
434 YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
435 YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
436 ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
437 ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
438
439 cout<<"XcoordMin "<<XcoordMin << " XcoordMax "<<XcoordMax <<endl;
440 cout<<"YcoordMin "<<YcoordMin << " YcoordMax "<<YcoordMax <<endl;
441 cout<<"ZcoordMin "<<ZcoordMin << " ZcoordMax "<<ZcoordMax <<endl;
442
443 /*double LxRVE, LyRVE, LzRVE;
444 LxRVE=XcoordMax-XcoordMin;
445 LyRVE=YcoordMax-YcoordMin;
446 LzRVE=ZcoordMax-ZcoordMin;*/
447
448 //Creating Prisims between triangles on -ve and +ve faces
449 typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
450 Tri_Hand_iterator Tri_Neg;
451 typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
452 xyzcoord_iterator Tri_Pos;
453 double XPos, YPos, ZPos;
454 //int count=0;
455
456 //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
457 int count=1;
458 int count1=0;
459 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
460
461// count1++;
462// cout<<"count1 = "<<count1<< endl;
463// if(count1 == 24) {
464
465
466 Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
467 //cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
468
469 //corresponding +ve triangle
470 if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
471 if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
472 if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
473 Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
474
475// cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
476// cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
477
478
479 //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
480 EntityHandle PrismNodes[6];
481 vector<EntityHandle> TriNodesNeg, TriNodesPos;
482 double CoordNodeNeg[9], CoordNodePos[9];
483 rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
484 rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
485 rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
486 rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
487 for(int ii=0; ii<3; ii++){
488 PrismNodes[ii]=TriNodesNeg[ii];
489 }
490
491// for(int ii=0; ii<3; ii++){
492// cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
493// }
494// for(int ii=0; ii<3; ii++){
495// cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
496// }
497
498 //Match exact nodes to each other to avoide the problem of twisted prisms
499 double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
500 for(int ii=0; ii<3; ii++){
501 if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
502 if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
503 if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
504 for(int jj=0; jj<3; jj++){
505 //Round nodal coordinates to 3 dicimal places only for comparison
506 XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
507 YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
508 ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
509
510 XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
511 XNodePos=round(XNodePos*roundfact)/roundfact;
512 YNodePos=round(YNodePos*roundfact)/roundfact;
513 ZNodePos=round(ZNodePos*roundfact)/roundfact;
514
515 if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
516 PrismNodes[3+ii]=TriNodesPos[jj];
517 break;
518 }
519 }
520 }
521 //prism nodes and their coordinates
522 double CoordNodesPrisms[18];
523 rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
524// for(int ii=0; ii<6; ii++){
525// cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
526// }
527// cout<<endl<<endl;
528
529 //insertion of individula prism element and its addition to range preriodic_prisms
530 EntityHandle PeriodicPrism;
531 rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
532 preriodic_prisms.insert(PeriodicPrism);
533
534// //to see individual prisms
535// Range Prism1;
536// Prism1.insert(PeriodicPrism);
537// EntityHandle out_meshset1;
538// rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
539// rval = moab.add_entities(out_meshset1,Prism1); CHKERRQ_MOAB(rval);
540// ostringstream sss;
541// sss << "Prism" << count << ".vtk"; count++;
542// rval = moab.write_file(sss.str().c_str(),"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
543//
544// string aaa;
545// cin>>aaa;
546//
547// }
548
549 }
550
551
552
553
554 //Adding the newely added prisms to the bit-level (which is known as seeding)
555 /*{
556 Range ents;
557 moab.get_adjacencies(preriodic_prisms,1,true,ents,moab::Interface::UNION);
558 moab.get_adjacencies(preriodic_prisms,2,true,ents,moab::Interface::UNION);
559 ierr = m_field.seed_ref_level(ents.subset_by_type(MBEDGE),problem_bit_level); CHKERRQ(ierr);
560 ierr = m_field.seed_ref_level(ents.subset_by_type(MBQUAD),problem_bit_level); CHKERRQ(ierr);
561 }*/
562 ierr = m_field.seed_ref_level(preriodic_prisms,problem_bit_level); CHKERRQ(ierr);
563 ierr = m_field.seed_finite_elements(preriodic_prisms); CHKERRQ(ierr);
564
565
566 //=======================================================================================================
567 //=======================================================================================================
568 //=======================================================================================================
569
570 EntityHandle out_meshset;
571 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
572 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
573
574 Range LatestRefinedTets;
575 rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
576
577 Range Interface_Periodic_Prisms;
578 rval = moab.get_entities_by_type(out_meshset, MBPRISM,Interface_Periodic_Prisms,true); CHKERRQ_MOAB(rval);
579 cout<<"Interface_Periodic_Prisms " <<Interface_Periodic_Prisms.size()<<endl;
580
581
582 //Fields
583 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
584 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
585 ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr); //For lagrange multipliers to control the periodic motion
586 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr); //To control the rigid body motion (3 Traslations and 3
587
588 //add entitities (by tets) to the field
589 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
590 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
591
592// ierr = PetscPrintf(PETSC_COMM_WORLD,"number of triangle in SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
593 ierr = m_field.add_ents_to_field_by_type(SurTrisNeg,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
594
595 //set app. order
596 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",bubble_order); CHKERRQ(ierr);
597 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
598 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
599 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
600
601 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
602 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
603 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
604 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
605
606 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
607 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
608 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
609
610 //build field
611 ierr = m_field.build_fields(); CHKERRQ(ierr);
612
613 //10 node tets
614 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
615 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material,0); CHKERRQ(ierr);
616
617 //FE
618 ierr = m_field.add_finite_element("PLASTIC"); CHKERRQ(ierr);
619 //Define rows/cols and element data
620 ierr = m_field.modify_finite_element_add_field_row("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
621 ierr = m_field.modify_finite_element_add_field_col("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
622 ierr = m_field.modify_finite_element_add_field_data("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
623 ierr = m_field.modify_finite_element_add_field_data("PLASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
624
625 Range newtets;
626 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
627 if(it->getName() != "MAT_PLASTIC") continue;
628 EntityHandle meshset = it->getMeshset();
629 int id = it->getMeshsetId();
630
631 rval = m_field.get_moab().get_entities_by_type(
632 meshset,MBTET,newtets,true
633 ); CHKERRQ_MOAB(rval);
634
635 cout<< "========================== newtets "<<newtets.size()<<endl;
636
637 //intersection of new and old tets for plastic
638 newtets = intersect(newtets,LatestRefinedTets);
639 ierr = m_field.seed_finite_elements(newtets); CHKERRQ(ierr);
640 }
641
642 ierr = m_field.add_ents_to_finite_element_by_type(newtets,MBTET,"PLASTIC"); CHKERRQ(ierr);
643 cout<< "========================== newtets "<<newtets.size()<<endl;
644
645
646 //================================================================================================================================
647 // Trans-Isotropic Yarns
648 //================================================================================================================================
649
650 NonlinearElasticElement trans_elastic(m_field,2);
651 trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
652 trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
653 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
654 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
655 bool trans_iso_blocks = false;
656 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
657 //Get block name
658 string name = it->getName();
659 if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
660 cout<<"================================ it is MAT_ELASTIC_TRANSISO "<<endl;
661 trans_iso_blocks = true;
662 int id = it->getMeshsetId();
663 Mat_Elastic_TransIso mydata;
664 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
665 tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
666 tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
667 //nu_p, nu_pz, E_p, E_z, G_zp
668 tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
669 tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
670 tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
671 tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
672 tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
673 tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
674 tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
675 tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
676 double shear_zp;
677 if(mydata.data.Shearzp!=0) {
678 shear_zp = mydata.data.Shearzp;
679 } else {
680 shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
681 }
682 tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
683 tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
684 //get tets from block where material is defined
685 EntityHandle meshset = it->getMeshset();
686 rval = m_field.get_moab().get_entities_by_type(
687 meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
688 ); CHKERRQ_MOAB(rval);
689
690// cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
691 //intersection of new and old tets for tran-iso part
692 trans_elastic.setOfBlocks[id].tEts = intersect(trans_elastic.setOfBlocks[id].tEts,LatestRefinedTets);
693// cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
694
695 //adding material to nonlinear class
696 trans_elastic.setOfBlocks[id].iD = id;
697 //note that material parameters are defined internally in material model
698 trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
699 trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
700 trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
701 trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
702 ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
703 }
704 }
705
706
707 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
708 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
709 ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
710 ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
711 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
712 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
713
714 if(m_field.check_field("MESH_NODE_POSITIONS")) {
715 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
716 }
717 for(
718 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
719 sit!=trans_elastic.setOfBlocks.end();sit++
720 ) {
721
722
723 cout<<" sit->second.tEts "<<sit->second.tEts.size()<<endl;
724 ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
725 }
726
727
728 //Rhs
729 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
730 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
731 if(m_field.check_field("MESH_NODE_POSITIONS")) {
732 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
733 }
734 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
735 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
736 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true));
737 trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpRhsPiolaKirchhoff("DISPLACEMENT",sit->second,trans_elastic.commonData));
738 }
739
740 //Lhs
741 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
742 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
743 if(m_field.check_field("MESH_NODE_POSITIONS")) {
744 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
745 }
746 sit = trans_elastic.setOfBlocks.begin();
747 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
748 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true));
749 trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx("DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData));
750 }
751
752
753
754
755 //================================================================================================================================
756 // Interface Cohesive elemetns
757 //================================================================================================================================
758 //FE Interface
759
760 ierr = m_field.add_finite_element("INTERFACE"); CHKERRQ(ierr);
761 ierr = m_field.modify_finite_element_add_field_row("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
762 ierr = m_field.modify_finite_element_add_field_col("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
763 ierr = m_field.modify_finite_element_add_field_data("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
764 ierr = m_field.modify_finite_element_add_field_data("INTERFACE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
765
766 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation> interface_materials;
767
768 //FIXME this in fact allow only for one type of interface,
769 //problem is Young Moduls in interface mayoung_modulusterial
770
771
772 double interface_beta, interface_ft, interface_Gf, interface_h;
773 ierr = PetscOptionsGetReal(0,"-interface_beta",&interface_beta,0); CHKERRQ(ierr);
774 ierr = PetscOptionsGetReal(0,"-interface_ft",&interface_ft,0); CHKERRQ(ierr);
775 ierr = PetscOptionsGetReal(0,"-interface_Gf",&interface_Gf,0); CHKERRQ(ierr);
776 ierr = PetscOptionsGetReal(0,"-interface_h",&interface_h,0); CHKERRQ(ierr);
777
778
779 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
780 cout << endl << *it << endl;
781 //Get block name
782 string name = it->getName();
783 if (name.compare(0,10,"MAT_INTERF") == 0) {
784 Mat_Interf mydata;
785 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
786 cout << mydata;
787 interface_materials.push_back(new CohesiveInterfaceElement::PhysicalEquation(m_field));
788// interface_materials.back().h = 0.134; // we used this because E0=35Gpa and Em=4.7Gpa
789 interface_materials.back().youngModulus = mydata.data.alpha;
790// interface_materials.back().beta = mydata.data.beta;
791// interface_materials.back().ft = mydata.data.ft;
792// interface_materials.back().Gf = mydata.data.Gf;
793
794 interface_materials.back().h = interface_h;
795 interface_materials.back().beta = interface_beta;
796 interface_materials.back().ft = interface_ft;
797 interface_materials.back().Gf = interface_Gf;
798
799 EntityHandle meshset = it->getMeshset();
800 Range tris; //All the triangles after splitting
801 rval = moab.get_entities_by_type(meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
802// cout<<"tris.size() "<< tris.size() <<endl;
803 Range ents3d;
804 rval = moab.get_adjacencies(tris,3,false,ents3d,moab::Interface::UNION); CHKERRQ_MOAB(rval);
805 interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
806// cout<<"unite( interface_prisms, ents3d.subset_by_type(MBPRISM) ) "<<unite( interface_prisms, ents3d.subset_by_type(MBPRISM) )<<endl;
807 }
808 }
809
810 ierr = m_field.seed_finite_elements(interface_materials.back().pRisms); CHKERRQ(ierr);
811// cout<<"interface_materials.back().pRisms.size() "<< interface_materials.back().pRisms <<endl;
812 ierr = m_field.add_ents_to_finite_element_by_type(interface_materials.back().pRisms,MBPRISM,"INTERFACE"); CHKERRQ(ierr);
813
814
815 // cout<<"young_modulus "<<young_modulus<<endl;
816 { //FIXME
817 boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit = interface_materials.begin();
818 for(; pit != interface_materials.end();pit++) {
819 pit->youngModulus = young_modulus; //Young's modulus of bulk elastic material (then for interface E0=E/h)
820 }
821 }
822
823 CohesiveInterfaceElement cohesive_elements(m_field); //make it golbal, so that its other uses can see it
824 ierr = cohesive_elements.addOps("DISPLACEMENT",interface_materials); CHKERRQ(ierr);
825 //================================================================================================================================
826
827
828 BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
829// BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
830
831 lagrangian_element_periodic.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",preriodic_prisms);
832
833
834 //Adding element LAGRANGE_ELEM_TRANS
835// lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
836 ierr = m_field.add_finite_element("LAGRANGE_ELEM_TRANS",MF_ZERO); CHKERRQ(ierr);
837 //============================================================================================================
838 //C row as Lagrange_mul_disp and col as DISPLACEMENT
839 ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
840 ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
841 //CT col as Lagrange_mul_disp and row as DISPLACEMENT
842 ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
843 ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
844 //data
845 ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
846 ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
847 //============================================================================================================
848
849
850 Range Tri103;
851 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
852 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
853
854 //total tri including old + new
855 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,Tri103,true); CHKERRQ_MOAB(rval);
856 ierr = PetscPrintf(PETSC_COMM_WORLD,"Old triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
857 //Only triangles in the new bit-level
858 Tri103 = intersect(AllTris, Tri103);
859 ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
860 ierr = m_field.add_ents_to_finite_element_by_type(Tri103,MBTRI,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
861 }
862 }
863
864// build finite elements
865 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
866 //build adjacencies
867 ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
868
869
870 DMType dm_name = "PLASTIC_PROB";
871 ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
872
873 DM dm;
874 ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
875 ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
876
877 //set dm datastruture which created mofem datastructures
878 ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,problem_bit_level); CHKERRQ(ierr);
879 ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
880 ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
881
882 //add elements to dm
883 ierr = DMMoFEMAddElement(dm,"PLASTIC"); CHKERRQ(ierr);
884 ierr = DMMoFEMAddElement(dm,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
885 ierr = DMMoFEMAddElement(dm,"INTERFACE"); CHKERRQ(ierr);
886 ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM"); CHKERRQ(ierr);
887 ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
888 ierr = DMSetUp(dm); CHKERRQ(ierr);
889
890 //create matrices
891 Vec F,D;
892 Mat Aij;
893 ierr = DMCreateGlobalVector_MoFEM(dm,&D); CHKERRQ(ierr);
894 ierr = VecDuplicate(D,&F); CHKERRQ(ierr);
895 ierr = DMCreateMatrix_MoFEM(dm,&Aij); CHKERRQ(ierr);
896 ierr = VecZeroEntries(D); CHKERRQ(ierr);
897 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
898 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
899 ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
900 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
901
902 vector<Vec> Fvec(6); //jthis will be used to caluclate the homogenised stiffness matix
903 for(int ii = 0;ii<6;ii++) {
904 ierr = VecDuplicate(D,&Fvec[ii]); CHKERRQ(ierr);
905 ierr = VecZeroEntries(Fvec[ii]); CHKERRQ(ierr);
906 ierr = VecGhostUpdateBegin(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
907 ierr = VecGhostUpdateEnd(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
908 }
909
910 //assemble Aij and F
911 SmallStrainPlasticity small_strain_plasticity(m_field);
912 {
913 PetscBool bbar = PETSC_TRUE;
914 ierr = PetscOptionsGetBool(0,"-my_bbar",&bbar,0); CHKERRQ(ierr);
915 small_strain_plasticity.commonData.bBar = bbar;
916 }
917 {
918 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
919 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
920 "DISPLACEMENT",small_strain_plasticity.commonData
921 )
922 );
923 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
924 new SmallStrainPlasticity::OpCalculateStress(
925 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
926 )
927 );
928 small_strain_plasticity.feRhs.getOpPtrVector().push_back(
929 new SmallStrainPlasticity::OpAssembleRhs(
930 "DISPLACEMENT",small_strain_plasticity.commonData
931 )
932 );
933 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
934 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
935 "DISPLACEMENT",small_strain_plasticity.commonData
936 )
937 );
938 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
939 new SmallStrainPlasticity::OpCalculateStress(
940 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
941 )
942 );
943 small_strain_plasticity.feLhs.getOpPtrVector().push_back(
944 new SmallStrainPlasticity::OpAssembleLhs(
945 "DISPLACEMENT",small_strain_plasticity.commonData
946 )
947 );
948 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
949 new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
950 "DISPLACEMENT",small_strain_plasticity.commonData
951 )
952 );
953 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
954 new SmallStrainPlasticity::OpCalculateStress(
955 m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
956 )
957 );
958 small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
959 new SmallStrainPlasticity::OpUpdate(
960 "DISPLACEMENT",small_strain_plasticity.commonData
961 )
962 );
963 ierr = small_strain_plasticity.createInternalVariablesTag(); CHKERRQ(ierr);
964 }
965
966
967
968 lagrangian_element_periodic.setRVEBCsOperatorsNonlinear("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,Fvec,F,given_strain);
969 BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
970 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element_periodic.setOfRVEBC);
971
972
973 TimeForceScale time_force_scale("-my_load_history",false);
974 lagrangian_element_periodic.methodsOp.push_back(new TimeForceScale("-my_macro_strian_history",false));
975
976 //Adding elements to DMSnes
977 //Rhs
978 ierr = DMMoFEMSNESSetFunction(dm,"PLASTIC",&small_strain_plasticity.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
979 ierr = DMMoFEMSNESSetFunction(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
980 ierr = DMMoFEMSNESSetFunction(dm,"INTERFACE",&cohesive_elements.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
981 ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
982 ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhsResidual,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
983
984 //Lhs
985 ierr = DMMoFEMSNESSetJacobian(dm,"PLASTIC",&small_strain_plasticity.feLhs,NULL,NULL); CHKERRQ(ierr);
986 ierr = DMMoFEMSNESSetJacobian(dm,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feLhs,NULL,NULL); CHKERRQ(ierr);
987 ierr = DMMoFEMSNESSetJacobian(dm,"INTERFACE",&cohesive_elements.feLhs,NULL,NULL); CHKERRQ(ierr);
988 ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
989 ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM_TRANS",&lagrangian_element_rigid_body_trans.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
990
991 // Create Time Stepping solver
992 SNES snes;
993 SnesCtx *snes_ctx;
994 SNES snes_one_step_only;
995
996
997 ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
998 //ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
999 ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
1000 ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
1001 ierr = SNESSetJacobian(snes,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
1002 ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
1003
1004 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_one_step_only); CHKERRQ(ierr);
1005 ierr = SNESSetFunction(snes_one_step_only,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
1006 ierr = SNESSetJacobian(snes_one_step_only,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
1007 ierr = SNESSetFromOptions(snes_one_step_only); CHKERRQ(ierr);
1008 double atol,rtol,stol;
1009 int maxit,maxf;
1010 SNESGetTolerances(snes_one_step_only,&atol,&rtol,&stol,&maxit,&maxf);
1011 maxit = 1;
1012 atol *= 0;
1013 SNESSetTolerances(snes_one_step_only,atol,rtol,stol,maxit,maxf);
1014 SNESLineSearch linesearch;
1015 SNESGetLineSearch(snes_one_step_only,&linesearch);
1016 SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);
1017
1018
1019
1020 PostProcVolumeOnRefinedMesh post_proc(m_field);
1021 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
1022 ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1023 ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1024 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1025
1026 // Volume calculation
1027 //==============================================================================================================================
1028 VolumeElementForcesAndSourcesCore vol_elem(m_field);
1029 Vec volume_vec;
1030 int volume_vec_ghost[] = { 0 };
1031 ierr = VecCreateGhost(
1032 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
1033 ); CHKERRQ(ierr);
1034 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
1035 vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1036
1037 ierr = m_field.loop_finite_elements("PLASTIC_PROB","PLASTIC", vol_elem); CHKERRQ(ierr);
1038 ierr = m_field.loop_finite_elements("PLASTIC_PROB","TRAN_ISOTROPIC_ELASTIC", vol_elem); CHKERRQ(ierr);
1039
1040 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
1041 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
1042 double rve_volume;
1043 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
1044
1045 // ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
1046 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
1047 ierr = VecDestroy(&volume_vec);
1048 //=============================================================================================================================
1049
1050 double final_time = 1,delta_time = 0.1;
1051 ierr = PetscOptionsGetReal(0,"-my_final_time",&final_time,0); CHKERRQ(ierr);
1052 ierr = PetscOptionsGetReal(0,"-my_delta_time",&delta_time,0); CHKERRQ(ierr);
1053 double delta_time0 = delta_time;
1054
1055// ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
1056// std::string wait;
1057// std::cin >> wait;
1058
1059 Vec D0;
1060 ierr = VecDuplicate(D,&D0); CHKERRQ(ierr);
1061
1062 int step = 0;
1063 double t = 0;
1064 SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
1065 for(;t<final_time;step++) {
1066 t += delta_time;
1067
1068// if(t>final_time){
1069// break;
1070// }
1071 PetscPrintf(
1072 PETSC_COMM_WORLD,"Step %d Time %6.4g final time %3.2g\n",step,t,final_time
1073 );
1074 //set time
1075 lagrangian_element_periodic.getLoopFeRVEBCRhs().ts_t = t;
1076 //solve problem at step
1077 ierr = VecAssemblyBegin(D);
1078 ierr = VecAssemblyEnd(D);
1079 ierr = VecCopy(D,D0); CHKERRQ(ierr);
1080 if(step == 0 || reason < 0) {
1081 ierr = SNESSetLagJacobian(snes,-2); CHKERRQ(ierr);
1082 } else {
1083 ierr = SNESSetLagJacobian(snes,-1); CHKERRQ(ierr);
1084 }
1085 ierr = SNESSetLagJacobian(snes_one_step_only,-2);
1086 ierr = SNESSolve(snes_one_step_only,PETSC_NULL,D); CHKERRQ(ierr);
1087 ierr = SNESGetConvergedReason(snes_one_step_only,&reason); CHKERRQ(ierr);
1088 if(reason<0) {
1089 ierr = SNESSolve(snes,PETSC_NULL,D); CHKERRQ(ierr);
1090 }
1091 int its;
1092 ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
1093
1094 ierr = PetscPrintf(
1095 PETSC_COMM_WORLD,"number of Newton iterations = %D\n",its
1096 ); CHKERRQ(ierr);
1097
1098 ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
1099
1100 if(reason<0) {
1101 t -= delta_time;
1102 delta_time *= 0.1;
1103 ierr = VecCopy(D0,D); CHKERRQ(ierr);
1104 } else {const int its_d = 60;
1105 const double gamma = 0.5;
1106 const double max_reudction = 1;
1107 const double min_reduction = 1e-1;
1108 double reduction;
1109 reduction = pow((double)its_d/(double)(its+1),gamma);
1110 if(delta_time >= max_reudction*delta_time0 && reduction > 1) {
1111 reduction = 1;
1112 } else if(delta_time <= min_reduction*delta_time0 && reduction < 1) {
1113 reduction = 1;
1114 }
1115
1116 ierr = PetscPrintf(
1117 PETSC_COMM_WORLD,
1118 "reduction delta_time = %6.4e delta_time = %6.4e\n",
1119 reduction,delta_time
1120 ); CHKERRQ(ierr);
1121 delta_time *= reduction;
1122 if(reduction>1 && delta_time < min_reduction*delta_time0) {
1123 delta_time = min_reduction*delta_time0;
1124 }
1125
1126 ierr = DMoFEMMeshToGlobalVector(
1127 dm,D,INSERT_VALUES,SCATTER_REVERSE
1128 ); CHKERRQ(ierr);
1129
1130 ierr = DMoFEMLoopFiniteElements(
1131 dm,"PLASTIC",&small_strain_plasticity.feUpdate
1132 ); CHKERRQ(ierr);
1133
1134
1135 ierr = DMoFEMLoopFiniteElements(
1136 dm,"INTERFACE",&cohesive_elements.feHistory
1137 ); CHKERRQ(ierr);
1138
1139 //Save data on mesh
1140 {
1141 ierr = DMoFEMLoopFiniteElements(
1142 dm,"PLASTIC",&post_proc
1143 ); CHKERRQ(ierr);
1144
1145 string out_file_name;
1146 {
1147 std::ostringstream stm;
1148 stm << "out_matrix_" << step << ".h5m";
1149 out_file_name = stm.str();
1150 }
1151
1152 ierr = PetscPrintf(
1153 PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
1154 ); CHKERRQ(ierr);
1155
1156 rval = post_proc.postProcMesh.write_file(
1157 out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
1158 ); CHKERRQ_MOAB(rval);
1159
1160
1161
1162 ierr = DMoFEMLoopFiniteElements(
1163 dm,"TRAN_ISOTROPIC_ELASTIC",&post_proc
1164 ); CHKERRQ(ierr);
1165 {
1166 std::ostringstream stm;
1167 stm << "out_fibres_" << step << ".h5m";
1168 out_file_name = stm.str();
1169 }
1170 ierr = PetscPrintf(
1171 PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
1172 ); CHKERRQ(ierr);
1173
1174 rval = post_proc.postProcMesh.write_file(
1175 out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
1176 ); CHKERRQ_MOAB(rval);
1177
1178 }
1179
1180 //===================================== Homgenised stress (for given strain) ====================================================
1181 VectorDouble StressHomo;
1182 StressHomo.resize(6);
1183 StressHomo.clear();
1184
1185 //calculate honmogenised stress
1186 //create a vector for 6 components of homogenized stress
1187 Vec stress_homo;
1188 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
1189 ierr = VecCreateGhost(
1190 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
1191 ); CHKERRQ(ierr);
1192
1193 lagrangian_element_periodic.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo);
1194
1195
1196 PetscScalar *avec;
1197 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
1198 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
1199 ierr = m_field.loop_finite_elements(
1200 "PLASTIC_PROB","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
1201 ); CHKERRQ(ierr);
1202 ierr = PetscOptionsGetReal(
1203 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
1204 ); CHKERRQ(ierr);
1205 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
1206 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
1207 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
1208 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1209 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1210
1211 for(int jj=0; jj<6; jj++) {
1212 StressHomo(jj) = avec[jj];
1213 }
1214
1215 double scale;
1216 ierr = time_force_scale.getForceScale(t,scale); CHKERRQ(ierr);
1217
1218 PetscPrintf(PETSC_COMM_WORLD,
1219 "Macro_Strain Homo_Stress Path "
1220 );
1221
1222 //cout<<"Macro-strain ";
1223 for(int ii=0; ii<6; ii++) {
1224 PetscPrintf(
1225 PETSC_COMM_WORLD,
1226 "%8.5e\t",
1227 t*given_strain(ii)
1228 );
1229 }
1230
1231 //cout<<"Homo stress ";
1232 for(int ii=0; ii<6; ii++) {
1233 PetscPrintf(
1234 PETSC_COMM_WORLD,
1235 "%8.5e\t",
1236 StressHomo(ii)
1237 );
1238 }
1239
1240 PetscPrintf(PETSC_COMM_WORLD,
1241 "\n"
1242 );
1243// //==============================================================================================================================
1244 }
1245 }
1246
1247 ierr = VecDestroy(&D0); CHKERRQ(ierr);
1248 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
1249 ierr = VecDestroy(&F); CHKERRQ(ierr);
1250 ierr = VecDestroy(&D); CHKERRQ(ierr);
1251
1252 }
1253 CATCH_ERRORS;
1254
1255
1256 MoFEM::Core::Finalize();
1257}
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
UNKNOWNNAME
@ UNKNOWNNAME
Definition definitions.h:171
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
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
MU
#define MU(E, NU)
Definition fem_tools.h:23
LAMBDA
#define LAMBDA(E, NU)
Definition fem_tools.h:22
F
@ F
Definition free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1123
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:497
MoFEM::DMMoFEMSNESSetFunction
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition DMMoFEM.cpp:718
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition DMMoFEM.cpp:114
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:523
MoFEM::DMMoFEMSNESSetJacobian
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition DMMoFEM.cpp:759
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1197
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:586
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition DMMoFEM.cpp:1167
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1094
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:535
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
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
out_file_name
char out_file_name[255]
Definition initial_diffusion.cpp:53
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::SnesMat
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition SnesCtx.cpp:139
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::SnesRhs
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition SnesCtx.cpp:27
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
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
young_modulus
double young_modulus
Young modulus.
Definition plastic.cpp:121
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:122
scale
double scale
Definition plastic.cpp:119
t
constexpr double t
plate stiffness
Definition plate.cpp:58
help
static char help[]
Definition rve_mech_plas_interface_periodic.cpp:71
BCs_RVELagrange_Periodic
Definition BCs_RVELagrange_Periodic.hpp:19
BCs_RVELagrange_Trac_Rigid_Trans
Definition BCs_RVELagrange_Trac_Rigid_Trans.hpp:18
Composite_xyzcoord
Definition homogenisation_periodic_atom_test.cpp:72
Face_CenPos_Handle
Definition homogenisation_periodic_atom_test.cpp:61
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
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::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_Interf
Linear interface data structure.
Definition MaterialBlocks.hpp:430
MoFEM::Mat_Interf::data
_data_ data
Definition MaterialBlocks.hpp:443
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:13
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition NonLinearElasticElement.hpp:557
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition HookeElement.hpp:27
PostProcVolumeOnRefinedMesh
Post processing.
Definition PostProcOnRefMesh.hpp:957
TimeForceScale
Force scale operator for reading two columns.
Definition TimeForceScale.hpp:18
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

Variable Documentation

◆ help

char help[]
static
Initial value:
=
"...\n"
"\n"

Definition at line 71 of file rve_mech_plas_interface_periodic.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow