v0.13.1
Loading...
Searching...
No Matches
rve_mechanical.cpp
Go to the documentation of this file.
1/** \file rve_mechanical.cpp
2 * \brief Calculates stiffness matrix for elastic RVE.
3
4 Three types of boundary conditions are implemented, i.e.
5 HOMOBCDISP, HOMOBCPERIODIC, HOMOBCTRAC, NITSCHE.
6
7 NITSHCE method allow to apply periodic boundary conditions
8 to arbitrary convex shape RVE.
9
10 */
11
12
13
14/* This file is part of MoFEM.
15 * MoFEM is free software: you can redistribute it and/or modify it under
16 * the terms of the GNU Lesser General Public License as published by the
17 * Free Software Foundation, either version 3 of the License, or (at your
18 * option) any later version.
19 *
20 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
21 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23 * License for more details.
24 *
25 * You should have received a copy of the GNU Lesser General Public
26 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
27
28#include <MoFEM.hpp>
29using namespace MoFEM;
30
31#include <Projection10NodeCoordsOnField.hpp>
32#include <petsctime.h>
33
34#include <boost/numeric/ublas/vector_proxy.hpp>
35
36#include <MethodForForceScaling.hpp>
37#include <TimeForceScale.hpp>
38
39#include <BCs_RVELagrange_Disp.hpp>
40#include <BCs_RVELagrange_Trac.hpp>
41#include <BCs_RVELagrange_Trac_Rigid_Rot.hpp>
42#include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
43#include <BCs_RVELagrange_Periodic.hpp>
44
45#ifndef WITH_ADOL_C
46 #error "MoFEM need to be compiled with ADOL-C"
47#endif
48
49#include <adolc/adolc.h>
50#include <NonLinearElasticElement.hpp>
51#include <Hooke.hpp>
52#include <boost/numeric/ublas/symmetric.hpp>
53#include <SmallTransverselyIsotropic.hpp>
54#include <VolumeCalculation.hpp>
55
56#include <boost/ptr_container/ptr_map.hpp>
57#include <PostProcOnRefMesh.hpp>
58#include <PostProcStresses.hpp>
59
60extern "C" {
61 void tetcircumcenter_tp(double a[3],double b[3],double c[3], double d[3],
62 double circumcenter[3],double *xi,double *eta,double *zeta);
63 void tricircumcenter3d_tp(double a[3],double b[3],double c[3],
64 double circumcenter[3],double *xi,double *eta);
65 //#include <triangle_ncc_rule.h>
66}
67
68// #include <NitscheMethod.hpp>
69#include <moab/AdaptiveKDTree.hpp>
70// #include <NitschePeriodicMethod.hpp>
71
72#include <algorithm>
73
74
75
76
77static char help[] = "...\n\n";
78
79enum HomoBCTypes {
80 HOMOBCDISP,
81 HOMOBCPERIODIC,
82 HOMOBCTRAC,
83 NITSCHE,
84 NOHOMOBC
85};
86
87const char *homo_bc_names[] = {
88 "disp",
89 "periodic",
90 "trac",
91 "nitsche",
92 "nohomobc"
93};
94
95//=================================================================================================================================
96//Define class and multindex container to store data for traiangles on the boundary of the RVE (it cannot be defined within main)
97//=================================================================================================================================
98
99struct Face_CenPos_Handle {
100 double xcoord, ycoord, zcoord;
101 const EntityHandle Tri_Hand;
102 Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand):xcoord(_xcoord),
103 ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
104};
105
106struct xcoord_tag {}; //tags to used in the multindex container
107struct ycoord_tag {};
108struct zcoord_tag {};
109struct Tri_Hand_tag {};
110struct Composite_xyzcoord {};
111
112typedef multi_index_container<
113Face_CenPos_Handle,
114indexed_by<
115ordered_non_unique<
116tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
117
118ordered_non_unique<
119tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
120
121ordered_non_unique<
122tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
123
124ordered_unique<
125tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
126
127ordered_unique<
128tag<Composite_xyzcoord>,
129composite_key<
130Face_CenPos_Handle,
131member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
132member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
133member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
134> > Face_CenPos_Handle_multiIndex;
135
136int main(int argc, char *argv[]) {
137
138 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
139
140 try {
141
142 moab::Core mb_instance;
143 moab::Interface& moab = mb_instance;
144 int rank;
145 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
146
147 //Reade parameters from line command
148 PetscBool flg = PETSC_TRUE;
149 char mesh_file_name[255];
150 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
151 if(flg != PETSC_TRUE) {
152 SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_FOUND,"*** ERROR -my_file (MESH FILE NEEDED)");
153 }
154
155 PetscInt order;
156 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
157 if(flg != PETSC_TRUE) {
158 order = 1;
159 }
160
161 PetscInt choise_value = NOHOMOBC;
162 ierr = PetscOptionsGetEList(
163 NULL,"-my_bc_type",homo_bc_names,NOHOMOBC,&choise_value,&flg
164 ); CHKERRQ(ierr);
165 if(flg != PETSC_TRUE) {
166 SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"*** Boundary conditions not set");
167 }
168
169 //Read mesh to MOAB
170 const char *option;
171 option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
172 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
173 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
174 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
175
176 //We need that for code profiling
177 PetscLogDouble t1,t2;
178 PetscLogDouble v1,v2;
179 ierr = PetscTime(&v1); CHKERRQ(ierr);
180 ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
181
182 //Create MoFEM (Joseph) database
183 MoFEM::Core core(moab);
184 MoFEM::Interface& m_field = core;
185
186 vector<BitRefLevel> bit_levels;
187 {
188 Tag th_meshset_info;
189 int def_meshset_info[2] = {0,0};
190 rval = moab.tag_get_handle(
191 "MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info
192 );
193 int meshset_data[2];
194 EntityHandle root = moab.get_root_set();
195 rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
196 if(meshset_data[0]==0) {
197 meshset_data[0] = 1;
198 rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
199
200 }
201 bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
202 }
203 BitRefLevel problem_bit_level = bit_levels.back();
204 ierr = m_field.seed_ref_level_3D(0,problem_bit_level); CHKERRQ(ierr);
205
206 // const clock_t begin_time = clock();
207 ierr = m_field.build_fields(); CHKERRQ(ierr);
208 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
209
210 Range preriodic_prisms;
211 if(choise_value == HOMOBCPERIODIC) {
212 //FIXME: Naming convention is not consistent in this section of code
213 //=======================================================================================================
214 //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
215 //=======================================================================================================
216 //Populating the Multi-index container with -ve triangles
217 Range SurTrisNeg;
218 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
219 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
220 Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
221 double TriCen[3], coords_Tri[9];
222
223
224 double roundfact=1e3;
225
226 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
227 const EntityHandle* conn_face; int num_nodes_Tri;
228
229 //get nodes attached to the face
230 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
231 //get nodal coordinates
232 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
233
234 //Find triangle centriod
235 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
236 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
237 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
238
239 //round values to roundfact disimal places
240 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
241 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
242 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
243
244 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
245 //fill the multi-index container with centriod coordinates and triangle handles
246 Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
247 // for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
248 // cout<<endl<<endl;
249 }
250
251 //Populating the Multi-index container with +ve triangles
252 Range SurTrisPos;
253 ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
254 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
255 for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end(); it++) {
256 const EntityHandle* conn_face; int num_nodes_Tri;
257
258 //get nodes attached to the face
259 rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
260 //get nodal coordinates
261 rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
262
263 //Find triangle centriod
264 TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
265 TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
266 TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
267
268 //round values to roundfact disimal places
269 TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
270 TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
271 TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
272 // cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << " TriCen[1]= "<< TriCen[1] << " TriCen[2]= "<< TriCen[2] <<endl;
273
274 //fill the multi-index container with centriod coordinates and triangle handles
275 Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
276 }
277
278 //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
279 double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
280 typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
281 typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
282 typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
283 Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
284 Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
285 Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
286
287 //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
288 XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin(); XcoordMin=XcoordMin_it->xcoord;
289 XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end(); XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
290 YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin(); YcoordMin=YcoordMin_it->ycoord;
291 YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end(); YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
292 ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin(); ZcoordMin=ZcoordMin_it->zcoord;
293 ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end(); ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
294
295 /*double LxRVE, LyRVE, LzRVE;
296 LxRVE=XcoordMax-XcoordMin;
297 LyRVE=YcoordMax-YcoordMin;
298 LzRVE=ZcoordMax-ZcoordMin;*/
299
300 //Creating Prisims between triangles on -ve and +ve faces
301 typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
302 Tri_Hand_iterator Tri_Neg;
303 typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
304 xyzcoord_iterator Tri_Pos;
305 double XPos, YPos, ZPos;
306 //int count=0;
307
308 //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
309 for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end(); it++) {
310
311 Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
312 // cout<<"xcoord= "<<Tri_iit->xcoord << " ycoord= "<< Tri_iit->ycoord << " ycoord= "<< Tri_iit->zcoord <<endl;
313
314 //corresponding +ve triangle
315 if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax; YPos=Tri_Neg->ycoord; ZPos=Tri_Neg->zcoord;};
316 if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=YcoordMax; ZPos=Tri_Neg->zcoord;};
317 if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord; YPos=Tri_Neg->ycoord; ZPos=ZcoordMax; };
318 Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
319 // cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << " Tri_Neg->ycoord "<< Tri_Neg->ycoord << " Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
320 // cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << " Tri_Pos->ycoord "<< Tri_Pos->ycoord << " Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
321
322
323 //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
324 EntityHandle PrismNodes[6];
325 vector<EntityHandle> TriNodesNeg, TriNodesPos;
326 double CoordNodeNeg[9], CoordNodePos[9];
327 rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
328 rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
329 rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg); MOAB_THROW(rval);
330 rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos); MOAB_THROW(rval);
331 for(int ii=0; ii<3; ii++){
332 PrismNodes[ii]=TriNodesNeg[ii];
333 }
334 // for(int ii=0; ii<3; ii++){
335 // cout<<"xcoord= "<<CoordNodeNeg[3*ii] << " ycoord= "<< CoordNodeNeg[3*ii+1] << " zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
336 // }
337 // for(int ii=0; ii<3; ii++){
338 // cout<<"xcoord= "<<CoordNodePos[3*ii] << " ycoord= "<< CoordNodePos[3*ii+1] << " zcoord= "<< CoordNodePos[3*ii+2] <<endl;
339 // }
340
341 //Match exact nodes to each other to avoide the problem of twisted prisms
342 double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
343 for(int ii=0; ii<3; ii++){
344 if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=CoordNodeNeg[3*ii+2];};
345 if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax; ZNodeNeg=CoordNodeNeg[3*ii+2];};
346 if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1]; ZNodeNeg=ZcoordMax;};
347 for(int jj=0; jj<3; jj++){
348 //Round nodal coordinates to 3 dicimal places only for comparison //round values to 3 disimal places
349 XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
350 YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
351 ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
352
353
354 XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
355 XNodePos=round(XNodePos*roundfact)/roundfact;
356 YNodePos=round(YNodePos*roundfact)/roundfact;
357 ZNodePos=round(ZNodePos*roundfact)/roundfact;
358
359 if(XNodeNeg==XNodePos && YNodeNeg==YNodePos && ZNodeNeg==ZNodePos){
360 PrismNodes[3+ii]=TriNodesPos[jj];
361 break;
362 }
363 }
364 }
365 //prism nodes and their coordinates
366 double CoordNodesPrisms[18];
367 rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms); MOAB_THROW(rval);
368 // for(int ii=0; ii<6; ii++){
369 // cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << " ycoord= "<< CoordNodesPrisms[3*ii+1] << " zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
370 // }
371 // cout<<endl<<endl;
372 //insertion of individula prism element and its addition to range preriodic_prisms
373 EntityHandle PeriodicPrism;
374 rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
375 preriodic_prisms.insert(PeriodicPrism);
376
377 }
378
379 //insertion of individual prism element and its addition to range preriodic_prisms
380 ierr = m_field.seed_ref_level(preriodic_prisms,problem_bit_level); CHKERRQ(ierr);
381
382 }
383
384 EntityHandle out_meshset;
385 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
386 // ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
387 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
388 Range LatestRefinedTets;
389 rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
390 Range LatestRefinedPrisms;
391 rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
392
393 Range prims_on_problem_bit_level;
394 ierr = m_field.get_entities_by_type_and_ref_level(
395 problem_bit_level,BitRefLevel().set(),MBPRISM,prims_on_problem_bit_level
396 ); CHKERRQ(ierr);
397 //to create meshset from range
398 EntityHandle meshset_prims_on_problem_bit_level;
399 rval = moab.create_meshset(MESHSET_SET,meshset_prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
400 rval = moab.add_entities(meshset_prims_on_problem_bit_level,prims_on_problem_bit_level); CHKERRQ_MOAB(rval);
401 ierr = m_field.seed_ref_level_MESHSET(meshset_prims_on_problem_bit_level,BitRefLevel().set()); CHKERRQ(ierr);
402
403 //Fields
404 int field_rank=3;
405 ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
406 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
407 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
408 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
409
410 //add entitities (by tets) to the field
411 if(choise_value == HOMOBCDISP) {
412 ierr = m_field.add_field("LAGRANGE_MUL_DISP",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
413 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
414 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
415 Range tris;
416 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
417 ierr = m_field.add_ents_to_field_by_type(tris,MBTRI,"LAGRANGE_MUL_DISP"); CHKERRQ(ierr);
418 }
419 }
420 }
421
422 if(choise_value == HOMOBCPERIODIC) {
423 ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
424 //Control 3 rigid body translations in x, y and z axis
425 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr);
426 // Setting up dummy no-field vertex
427 EntityHandle no_field_vertex;
428 {
429 const double coords[] = {0,0,0};
430 rval = m_field.get_moab().create_vertex(coords,no_field_vertex); CHKERRQ_MOAB(rval);
431 Range range_no_field_vertex;
432 range_no_field_vertex.insert(no_field_vertex);
433 ierr = m_field.seed_ref_level(range_no_field_vertex,BitRefLevel().set()); CHKERRQ(ierr);
434 EntityHandle meshset;
435 meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_TRANS");
436 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
437 }
438 Range surface_negative;
439 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,surface_negative,true); CHKERRQ(ierr);
440 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",surface_negative.size()); CHKERRQ(ierr);
441 ierr = m_field.add_ents_to_field_by_type(surface_negative,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
442 }
443
444 if(choise_value == HOMOBCTRAC) {
445 ierr = m_field.add_field("LAGRANGE_MUL_TRAC",NOFIELD,NOBASE,6); CHKERRQ(ierr);
446 //Control 3 rigid body translations in x, y and z axis
447 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr);
448 //Controla 3 rigid body rotations about x, y and z axis
449 ierr = m_field.add_field("LAGRANGE_MUL_RIGID_ROT",NOFIELD,NOBASE,3); CHKERRQ(ierr);
450 EntityHandle no_field_vertex;
451 {
452 const double coords[] = {0,0,0};
453 rval = m_field.get_moab().create_vertex(coords,no_field_vertex); CHKERRQ_MOAB(rval);
454 Range range_no_field_vertex;
455 range_no_field_vertex.insert(no_field_vertex);
456 ierr = m_field.seed_ref_level(range_no_field_vertex,BitRefLevel().set()); CHKERRQ(ierr);
457 EntityHandle meshset;
458 meshset = m_field.get_field_meshset("LAGRANGE_MUL_TRAC");
459 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
460 meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_TRANS");
461 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
462 meshset = m_field.get_field_meshset("LAGRANGE_MUL_RIGID_ROT");
463 rval = m_field.get_moab().add_entities(meshset,range_no_field_vertex); CHKERRQ_MOAB(rval);
464 }
465 }
466
467 //set app. order
468 //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
469 ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
470 ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
471 ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
472 ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
473
474 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
475 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
476 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
477 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
478
479 PetscBool fo_boundary = PETSC_FALSE;
480 ierr = PetscOptionsGetBool(PETSC_NULL,"-my_fo_boundary",&fo_boundary,PETSC_NULL); CHKERRQ(ierr);
481 if(fo_boundary) {
482 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
483 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
484 Range tris;
485 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
486 Range tris_edges;
487 rval = moab.get_adjacencies(tris,1,false,tris_edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
488 ierr = m_field.set_field_order(tris,"DISPLACEMENT",1); CHKERRQ(ierr);
489 ierr = m_field.set_field_order(tris_edges,"DISPLACEMENT",1); CHKERRQ(ierr);
490 }
491 }
492 }
493
494 if(choise_value == HOMOBCDISP) {
495 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
496 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
497 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_DISP",1); CHKERRQ(ierr);
498 }
499
500 if(choise_value == HOMOBCPERIODIC) {
501 ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
502 ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
503 ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
504 }
505
506 //build field
507 ierr = m_field.build_fields(); CHKERRQ(ierr);
508
509 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
510 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
511
512 boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>());
513 boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>());
514
515 NonlinearElasticElement iso_elastic(m_field,1);
516 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|MAT_ELASTICSET,it)) {
517 if(it->getName() != "MAT_ELASTIC_1") continue;
518 Mat_Elastic mydata;
519 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
520 int id = it->getMeshsetId();
521 EntityHandle meshset = it->getMeshset();
522 rval = m_field.get_moab().get_entities_by_type(
523 meshset,MBTET,iso_elastic.setOfBlocks[id].tEts,true
524 ); CHKERRQ_MOAB(rval);
525 iso_elastic.setOfBlocks[id].iD = id;
526 iso_elastic.setOfBlocks[id].E = mydata.data.Young;
527 iso_elastic.setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
528 iso_elastic.setOfBlocks[id].materialDoublePtr = hooke_double;
529 iso_elastic.setOfBlocks[id].materialAdoublePtr = hooke_adouble;
530 ierr = m_field.seed_finite_elements(iso_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
531 }
532 ierr = iso_elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
533 ierr = iso_elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
534 if(m_field.check_field("POTENTIAL_FIELD")) {
535 ierr = m_field.modify_finite_element_add_field_data("ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
536 }
537
538 NonlinearElasticElement trans_elastic(m_field,2);
539 trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
540 trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
541 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
542 std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
543 bool trans_iso_blocks = false;
544 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
545 //Get block name
546 string name = it->getName();
547 if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
548 trans_iso_blocks = true;
549 int id = it->getMeshsetId();
550 Mat_Elastic_TransIso mydata;
551 ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
552 tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
553 tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
554 //nu_p, nu_pz, E_p, E_z, G_zp
555 tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
556 tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
557 tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
558 tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
559 tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
560 tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
561 tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
562 tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
563 double shear_zp;
564 if(mydata.data.Shearzp!=0) {
565 shear_zp = mydata.data.Shearzp;
566 } else {
567 shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
568 }
569 tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
570 tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
571 //get tets from block where material is defined
572 EntityHandle meshset = it->getMeshset();
573 rval = m_field.get_moab().get_entities_by_type(
574 meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
575 ); CHKERRQ_MOAB(rval);
576 //adding material to nonlinear class
577 trans_elastic.setOfBlocks[id].iD = id;
578 //note that material parameters are defined internally in material model
579 trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
580 trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
581 trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
582 trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
583 ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
584 }
585 }
586 if(trans_iso_blocks) {
587 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
588 ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
589 ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
590 ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
591 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
592 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
593 if(m_field.check_field("MESH_NODE_POSITIONS")) {
594 ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
595 }
596 for(
597 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
598 sit!=trans_elastic.setOfBlocks.end();sit++
599 ) {
600 ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
601 }
602 }
603 if(trans_iso_blocks) {
604 //Rhs
605 trans_elastic.feRhs.getOpPtrVector().push_back(
606 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData)
607 );
608 trans_elastic.feRhs.getOpPtrVector().push_back(
609 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData)
610 );
611 if(m_field.check_field("MESH_NODE_POSITIONS")) {
612 trans_elastic.feRhs.getOpPtrVector().push_back(
613 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData)
614 );
615 }
616 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
617 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
618 trans_elastic.feRhs.getOpPtrVector().push_back(
619 new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
620 "DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true
621 )
622 );
623 trans_elastic.feRhs.getOpPtrVector().push_back(
624 new NonlinearElasticElement::OpRhsPiolaKirchhoff(
625 "DISPLACEMENT",sit->second,trans_elastic.commonData
626 )
627 );
628 }
629
630 //Lhs
631 trans_elastic.feLhs.getOpPtrVector().push_back(
632 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData)
633 );
634 trans_elastic.feLhs.getOpPtrVector().push_back(
635 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData)
636 );
637 if(m_field.check_field("MESH_NODE_POSITIONS")) {
638 trans_elastic.feLhs.getOpPtrVector().push_back(
639 new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData)
640 );
641 }
642 sit = trans_elastic.setOfBlocks.begin();
643 for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
644 trans_elastic.feLhs.getOpPtrVector().push_back(
645 new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
646 "DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true
647 )
648 );
649 trans_elastic.feLhs.getOpPtrVector().push_back(
650 new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(
651 "DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData
652 )
653 );
654 }
655 }
656
657 BCs_RVELagrange_Disp lagrangian_element_disp(m_field);
658 if(choise_value == HOMOBCDISP) {
659 lagrangian_element_disp.addLagrangiangElement(
660 "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS"
661 );
662 }
663
664 BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
665 BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
666 BCs_RVELagrange_Trac_Rigid_Rot lagrangian_element_rigid_body_rot(m_field);
667 if(choise_value == HOMOBCTRAC) {
668 lagrangian_element_trac.addLagrangiangElement(
669 "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS"
670 );
671 lagrangian_element_trac.addLagrangiangElement(
672 "LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS"
673 );
674 lagrangian_element_trac.addLagrangiangElement(
675 "LAGRANGE_ELEM_ROT","DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT","MESH_NODE_POSITIONS"
676 );
677 }
678
679 BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
680 if(choise_value == HOMOBCPERIODIC) {
681 lagrangian_element_periodic.addLagrangiangElement(
682 "LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",preriodic_prisms
683 );
684 lagrangian_element_trac.addLagrangiangElement(
685 "LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS"
686 );
687 }
688
689 struct MinMaxNodes {
690 enum MINAMX { C0,MAXLAST };
691 EntityHandle entMinMax[MAXLAST];
692 ublas::vector<int> rowIndices;
693 VectorDouble rHs[6];
694 MinMaxNodes() {
695 rowIndices.resize(3*MAXLAST);
696 for(int ii = 0;ii<6;ii++) {
697 rHs[ii].resize(3*MAXLAST);
698 }
699 }
700 };
701 MinMaxNodes minMaxNodes;
702
703 if(choise_value == NITSCHE) { // Condensed traction/periodc BC
704 ierr = m_field.add_finite_element("SURFACE_ELEMENTS"); CHKERRQ(ierr);
705 ierr = m_field.modify_finite_element_add_field_row("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
706 ierr = m_field.modify_finite_element_add_field_col("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
707 ierr = m_field.modify_finite_element_add_field_data("SURFACE_ELEMENTS","DISPLACEMENT"); CHKERRQ(ierr);
708 ierr = m_field.modify_finite_element_add_field_data("SURFACE_ELEMENTS","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
709 EntityHandle condensed_traction_element_meshset;
710 rval = moab.create_meshset(MESHSET_TRACK_OWNER,condensed_traction_element_meshset); CHKERRQ(ierr);
711 Range nodes;
712 for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
713 if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
714 Range tris;
715 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
716 ierr = m_field.add_ents_to_finite_element_by_type(tris,MBTRI,"SURFACE_ELEMENTS"); CHKERRQ(ierr);
717 Range tris_nodes;
718 rval = moab.get_connectivity(tris,nodes,true); CHKERRQ_MOAB(rval);
719 nodes.merge(tris_nodes);
720 }
721 }
722
723 {
724 VectorDouble x,y,z;
725 x.resize(nodes.size(),false);
726 y.resize(nodes.size(),false);
727 z.resize(nodes.size(),false);
728 rval = moab.get_coords(nodes,&x[0],&y[0],&z[0]); CHKERRQ_MOAB(rval);
729 int bc_nb = 0;
730 for(int sx = -1; sx<=+1; sx+=2) {
731 for(int sy = -1; sy<=+1; sy+=2) {
732 for(int sz = -1; sz<=+1; sz+=2) {
733 if(bc_nb == MinMaxNodes::MAXLAST) break;
734 VectorDouble dist_up_right;
735 dist_up_right.resize(x.size(),false);
736 dist_up_right.clear();
737 for(unsigned int nn = 0;nn<x.size();nn++) {
738 if(
739 ((sx*x[nn])>0)&&
740 ((sy*y[nn])>0)&&
741 ((sz*z[nn])>0)
742 ) {
743 dist_up_right[nn] = sx*x[nn]+sy*y[nn]+sz*z[nn];
744 }
745 }
746 VectorDouble::iterator dist_it;
747 dist_it = max_element(dist_up_right.begin(),dist_up_right.end());
748 minMaxNodes.entMinMax[bc_nb++] = nodes[std::distance(dist_up_right.begin(),dist_it)];
749 }
750 }
751 }
752 }
753
754 }
755
756 //build finite elements
757 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
758 //build adjacencies
759 ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
760
761 //define problems
762 ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
763 //set finite elements for problem
764 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
765 if(trans_iso_blocks) {
766 ierr = m_field.modify_problem_add_finite_element(
767 "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC"
768 ); CHKERRQ(ierr);
769 }
770 if(choise_value == HOMOBCDISP) {
771 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
772 }
773 if(choise_value == HOMOBCTRAC) {
774 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
775 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
776 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT"); CHKERRQ(ierr);
777 }
778 if(choise_value == HOMOBCPERIODIC) {
779 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
780 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
781 }
782 if(choise_value == NITSCHE) {
783 ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","SURFACE_ELEMENTS"); CHKERRQ(ierr);
784 }
785
786 //set refinement level for problem
787 ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",problem_bit_level); CHKERRQ(ierr);
788 //build problem
789 ierr = m_field.build_problems(); CHKERRQ(ierr);
790
791 /****/
792 //mesh partitioning
793
794 //partition
795 ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
796 ierr = m_field.partition_finite_elements(
797 "ELASTIC_MECHANICS",false,0,m_field.get_comm_size() // build elements on all procs
798 ); CHKERRQ(ierr);
799 //what are ghost nodes, see Petsc Manual
800 ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
801
802 //create matrices
803 Vec D;
804 vector<Vec> F(6);
805 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
806 for(int ii = 1;ii<6;ii++) {
807 ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
808 }
809 ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
810
811 Mat Aij;
812 ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
813 ierr = MatSetOption(Aij,MAT_STRUCTURALLY_SYMMETRIC,PETSC_TRUE); CHKERRQ(ierr);
814 ierr = MatSetOption(Aij,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE); CHKERRQ(ierr);
815 ierr = MatSetOption(Aij,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE); CHKERRQ(ierr);
816 ierr = MatSetOption(Aij,MAT_USE_INODES,PETSC_TRUE); CHKERRQ(ierr);
817 ierr = MatSetOption(Aij,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE); CHKERRQ(ierr);
818 ierr = MatSetOption(Aij,MAT_KEEP_NONZERO_PATTERN,PETSC_FALSE); CHKERRQ(ierr);
819
820 /*{
821 ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
822 std::string wait;
823 std::cin >> wait;
824 }*/
825
826 ierr = VecZeroEntries(D); CHKERRQ(ierr);
827 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
828 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
829 ierr = m_field.set_global_ghost_vector(
830 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
831 ); CHKERRQ(ierr);
832 for(int ii = 0;ii<6;ii++) {
833 ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
834 ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
835 ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
836 }
837 ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
838
839 // NitscheMethod::BlockData nitsche_block_data;
840 // NitscheMethod::CommonData nitsche_common_data;
841 // PeriodicNitscheConstrains::CommonData periodic_common_data;
842 // PeriodicNitscheConstrains::MyNitscheVolume nitsche_element_iso(
843 // m_field,nitsche_block_data,nitsche_common_data,periodic_common_data
844 // );
845 // PeriodicNitscheConstrains::MyNitscheVolume nitsche_element_trans(
846 // m_field,nitsche_block_data,nitsche_common_data,periodic_common_data
847 // );
848
849 // if(choise_value == NITSCHE) {
850
851 // nitsche_block_data.faceElemName = "SURFACE_ELEMENTS";
852 // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
853 // if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
854 // Range tris;
855 // rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
856 // nitsche_block_data.fAces.merge(tris);
857 // }
858 // }
859
860 // for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
861 // if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
862 // Range tris;
863 // rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
864 // periodic_common_data.skinFaces.merge(tris);
865 // }
866 // }
867
868 // nitsche_block_data.gamma = 1e-7;
869 // nitsche_block_data.phi = -1;
870 // periodic_common_data.eps = 0;
871 // ierr = PetscOptionsGetReal(
872 // PETSC_NULL,"-my_gamma",&nitsche_block_data.gamma,PETSC_NULL
873 // ); CHKERRQ(ierr);
874 // ierr = PetscOptionsGetReal(
875 // PETSC_NULL,"-my_phi",&nitsche_block_data.phi,PETSC_NULL
876 // ); CHKERRQ(ierr);
877 // ierr = PetscOptionsGetReal(
878 // PETSC_NULL,"-my_eps",&periodic_common_data.eps,PETSC_NULL
879 // ); CHKERRQ(ierr);
880 // ierr = PetscPrintf(
881 // PETSC_COMM_WORLD,
882 // "Nitsche method gamma = %4.2e phi = %2.1f eps = %4.2e\n",
883 // nitsche_block_data.gamma,nitsche_block_data.phi,periodic_common_data.eps
884 // ); CHKERRQ(ierr);
885
886 // for(
887 // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
888 // mit!=iso_elastic.setOfBlocks.end();
889 // mit++
890 // ) {
891 // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
892 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
893 // nitsche_element_iso.snes_B = Aij;
894
895 // nitsche_element_iso.getOpPtrVector().push_back(
896 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
897 // "DISPLACEMENT",elastic_common_data
898 // )
899 // );
900 // nitsche_element_iso.getOpPtrVector().push_back(
901 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
902 // "MESH_NODE_POSITIONS",elastic_common_data
903 // )
904 // );
905 // nitsche_element_iso.getOpPtrVector().push_back(
906 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
907 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,true,false,true
908 // )
909 // );
910 // nitsche_element_iso.getOpPtrVector().push_back(
911 // new PeriodicNitscheConstrains::OpLhsPeriodicNormal(
912 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
913 // elastic_block_data,elastic_common_data,
914 // periodic_common_data
915 // )
916 // );
917 // nitsche_element_iso.getOpPtrVector().push_back(
918 // new PeriodicNitscheConstrains::OpRhsPeriodicNormal(
919 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
920 // elastic_block_data,elastic_common_data,
921 // periodic_common_data,
922 // F
923 // )
924 // );
925 // nitsche_element_iso.getOpPtrVector().push_back(
926 // new NitscheMethod::OpLhsNormal(
927 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
928 // elastic_block_data,elastic_common_data,true
929 // )
930 // );
931
932 // // this is to get data on opposite element
933 // nitsche_element_iso.periodicVolume.getOpPtrVector().clear();
934 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
935 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
936 // "DISPLACEMENT",elastic_common_data
937 // )
938 // );
939 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
940 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
941 // "MESH_NODE_POSITIONS",elastic_common_data
942 // )
943 // );
944 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
945 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
946 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,true,false,true
947 // )
948 // );
949 // nitsche_element_iso.periodicVolume.getOpPtrVector().push_back(
950 // new PeriodicNitscheConstrains::OpGetVolumeData(
951 // elastic_common_data,
952 // periodic_common_data
953 // )
954 // );
955 // periodic_common_data.volumeElemName = "ELASTIC";
956
957 // nitsche_element_iso.addToRule = 1;
958 // ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",nitsche_element_iso); CHKERRQ(ierr);
959 // }
960
961 // for(
962 // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
963 // mit!=trans_elastic.setOfBlocks.end();
964 // mit++
965 // ) {
966 // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
967 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
968 // nitsche_element_trans.snes_B = Aij;
969
970 // nitsche_element_trans.getOpPtrVector().push_back(
971 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
972 // "DISPLACEMENT",elastic_common_data
973 // )
974 // );
975 // nitsche_element_trans.getOpPtrVector().push_back(
976 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
977 // "MESH_NODE_POSITIONS",elastic_common_data
978 // )
979 // );
980 // nitsche_element_trans.getOpPtrVector().push_back(
981 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
982 // "POTENTIAL_FIELD",elastic_common_data
983 // )
984 // );
985 // nitsche_element_trans.getOpPtrVector().push_back(
986 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
987 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,true,false,true
988 // )
989 // );
990 // nitsche_element_trans.getOpPtrVector().push_back(
991 // new PeriodicNitscheConstrains::OpLhsPeriodicNormal(
992 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
993 // elastic_block_data,elastic_common_data,
994 // periodic_common_data
995 // )
996 // );
997 // nitsche_element_trans.getOpPtrVector().push_back(
998 // new PeriodicNitscheConstrains::OpRhsPeriodicNormal(
999 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1000 // elastic_block_data,elastic_common_data,
1001 // periodic_common_data,
1002 // F
1003 // )
1004 // );
1005 // nitsche_element_trans.getOpPtrVector().push_back(
1006 // new NitscheMethod::OpLhsNormal(
1007 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1008 // elastic_block_data,elastic_common_data,true
1009 // )
1010 // );
1011
1012 // // this is to get data on opposite element
1013 // nitsche_element_trans.periodicVolume.getOpPtrVector().clear();
1014 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1015 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1016 // "DISPLACEMENT",elastic_common_data
1017 // )
1018 // );
1019 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1020 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1021 // "POTENTIAL_FIELD",elastic_common_data
1022 // )
1023 // );
1024 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1025 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1026 // "MESH_NODE_POSITIONS",elastic_common_data
1027 // )
1028 // );
1029 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1030 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1031 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,true,false,true
1032 // )
1033 // );
1034 // nitsche_element_trans.periodicVolume.getOpPtrVector().push_back(
1035 // new PeriodicNitscheConstrains::OpGetVolumeData(
1036 // elastic_common_data,
1037 // periodic_common_data
1038 // )
1039 // );
1040 // periodic_common_data.volumeElemName = "TRAN_ISOTROPIC_ELASTIC";
1041
1042 // nitsche_element_trans.addToRule = 1;
1043 // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1044 // ierr = m_field.loop_finite_elements(
1045 // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",nitsche_element_trans
1046 // ); CHKERRQ(ierr);
1047 // }
1048
1049
1050 // }
1051
1052 Vec volume_vec;
1053 int volume_vec_ghost[] = { 0 };
1054 ierr = VecCreateGhost(
1055 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
1056 ); CHKERRQ(ierr);
1057 ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
1058
1059 iso_elastic.getLoopFeLhs().getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1060 trans_elastic.getLoopFeLhs().getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
1061
1062 //iso_elastic element matrix
1063 iso_elastic.getLoopFeLhs().snes_x = D;
1064 iso_elastic.getLoopFeLhs().snes_B = Aij;
1065 trans_elastic.getLoopFeLhs().snes_x = D;
1066 trans_elastic.getLoopFeLhs().snes_B = Aij;
1067 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",iso_elastic.getLoopFeLhs()); CHKERRQ(ierr);
1068 if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1069 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",trans_elastic.getLoopFeLhs()); CHKERRQ(ierr);
1070
1071 if(choise_value == HOMOBCDISP) {
1072 lagrangian_element_disp.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",Aij,F);
1073 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
1074 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
1075 }
1076 if(choise_value == HOMOBCTRAC) {
1077 lagrangian_element_trac.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",Aij,F);
1078 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
1079 ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
1080 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators(
1081 "DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_trac.setOfRVEBC
1082 );
1083 ierr = m_field.loop_finite_elements(
1084 "ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()
1085 ); CHKERRQ(ierr);
1086 lagrangian_element_rigid_body_rot.setRVEBCsRigidBodyRotOperators(
1087 "DISPLACEMENT","LAGRANGE_MUL_RIGID_ROT",Aij,lagrangian_element_trac.setOfRVEBC
1088 );
1089 ierr = m_field.loop_finite_elements(
1090 "ELASTIC_MECHANICS","LAGRANGE_ELEM_ROT",lagrangian_element_rigid_body_rot.getLoopFeRVEBCLhs()
1091 ); CHKERRQ(ierr);
1092 }
1093 if(choise_value == HOMOBCPERIODIC) {
1094 lagrangian_element_periodic.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,F);
1095 ierr = m_field.loop_finite_elements(
1096 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCLhs()
1097 ); CHKERRQ(ierr);
1098 ierr = m_field.loop_finite_elements(
1099 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCRhs()
1100 ); CHKERRQ(ierr);
1101 lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators(
1102 "DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij,lagrangian_element_periodic.setOfRVEBC
1103 );
1104 ierr = m_field.loop_finite_elements(
1105 "ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()
1106 ); CHKERRQ(ierr);
1107 }
1108
1109 for(int ii = 0;ii<6;ii++) {
1110 ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
1111 ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
1112 ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
1113 ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
1114 }
1115 ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1116 ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
1117
1118 {
1119 //ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
1120 //std::string wait;
1121 //std::cin >> wait;
1122 }
1123
1124 // if(choise_value == NITSCHE) {
1125 // if(periodic_common_data.eps==0) {
1126 // const Problem *problem_ptr;
1127 // ierr = m_field.get_problem("ELASTIC_MECHANICS",&problem_ptr); CHKERRQ(ierr);
1128 // for(int nn = 0;nn!=MinMaxNodes::MAXLAST;nn++) {
1129 // for(
1130 // _IT_NUMEREDDOF_ROW_BY_ENT_FOR_LOOP_(
1131 // problem_ptr,minMaxNodes.entMinMax[nn],dof
1132 // )
1133 // ) {
1134 // minMaxNodes.rowIndices[3*nn+dof->get()->getDofCoeffIdx()]
1135 // = dof->get()->getPetscGlobalDofIdx();
1136 // }
1137 // }
1138 // VectorDouble coords;
1139 // int nb_bcs = 3*MinMaxNodes::MAXLAST;
1140 // coords.resize(nb_bcs,false);
1141 // rval = moab.get_coords(&minMaxNodes.entMinMax[0],MinMaxNodes::MAXLAST,&coords[0]);
1142 // VectorDouble strain;
1143 // strain.resize(6,false);
1144 // MatrixDouble mat_d;
1145 // for(int ii = 0;ii<6;ii++) {
1146 // minMaxNodes.rHs[ii].clear();
1147 // strain.clear();
1148 // strain[ii] = 1;
1149 // for(int nn = 0;nn<MinMaxNodes::MAXLAST;nn++) {
1150 // PeriodicNitscheConstrains::OpRhsPeriodicNormal::calcualteDMatrix(
1151 // VectorAdaptor(3,ublas::shallow_array_adaptor<double>(3,&coords[3*nn])),
1152 // mat_d
1153 // );
1154 // VectorAdaptor rhs(3,ublas::shallow_array_adaptor<double>(3,&minMaxNodes.rHs[ii][3*nn]));
1155 // noalias(rhs) = prod(mat_d,strain);
1156 // }
1157 // }
1158 // for(int ii = 0;ii<6;ii++) {
1159 // ierr = VecSetValues(
1160 // F[ii],nb_bcs,&minMaxNodes.rowIndices[0],&minMaxNodes.rHs[ii][0],INSERT_VALUES
1161 // ); CHKERRQ(ierr);
1162 // ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
1163 // ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
1164 // }
1165 // ierr = MatZeroRows(
1166 // Aij,nb_bcs,&minMaxNodes.rowIndices[0],1,PETSC_NULL,PETSC_NULL
1167 // ); CHKERRQ(ierr);
1168 // }
1169 // }
1170
1171 ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
1172 ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
1173 double rve_volume;
1174 ierr = VecSum(volume_vec,&rve_volume); CHKERRQ(ierr);
1175 ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
1176
1177 ierr = VecDestroy(&volume_vec);
1178
1179 // Solver
1180 KSP solver;
1181 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
1182 ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
1183 ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
1184 ierr = KSPSetUp(solver); CHKERRQ(ierr);
1185
1186 MatrixDouble Dmat;
1187 Dmat.resize(6,6);
1188 Dmat.clear();
1189
1190 //create a vector for 6 components of homogenized stress
1191 Vec stress_homo;
1192 int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
1193 NonlinearElasticElement::MyVolumeFE ave_stress_iso(m_field);
1194 NonlinearElasticElement::MyVolumeFE ave_stress_trans(m_field);
1195 PetscBool stress_by_boundary_integral = PETSC_FALSE;
1196 ierr = VecCreateGhost(
1197 PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
1198 ); CHKERRQ(ierr);
1199 switch(choise_value) {
1200 case HOMOBCDISP:
1201 lagrangian_element_disp.setRVEBCsHomoStressOperators(
1202 "DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",stress_homo
1203 );
1204 break;
1205 case HOMOBCTRAC:
1206 lagrangian_element_trac.setRVEBCsHomoStressOperators(
1207 "DISPLACEMENT","LAGRANGE_MUL_TRAC","MESH_NODE_POSITIONS",stress_homo
1208 );
1209 break;
1210 case HOMOBCPERIODIC:
1211 lagrangian_element_periodic.setRVEBCsHomoStressOperators(
1212 "DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo
1213 );
1214 break;
1215 case NITSCHE:
1216 // if(stress_by_boundary_integral) {
1217 // nitsche_element_iso.getOpPtrVector().clear();
1218 // nitsche_element_trans.getOpPtrVector().clear();
1219 // nitsche_element_iso.periodicVolume.getOpPtrVector().clear();
1220 // nitsche_element_trans.periodicVolume.getOpPtrVector().clear();
1221 // for(
1222 // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
1223 // mit!=iso_elastic.setOfBlocks.end();
1224 // mit++
1225 // ) {
1226 // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
1227 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1228 // nitsche_element_iso.getOpPtrVector().push_back(
1229 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1230 // "DISPLACEMENT",elastic_common_data
1231 // )
1232 // );
1233 // nitsche_element_iso.getOpPtrVector().push_back(
1234 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1235 // "MESH_NODE_POSITIONS",elastic_common_data
1236 // )
1237 // );
1238 // nitsche_element_iso.getOpPtrVector().push_back(
1239 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1240 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,false,false,true
1241 // )
1242 // );
1243 // nitsche_element_iso.getOpPtrVector().push_back(
1244 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral(
1245 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1246 // elastic_block_data,elastic_common_data,stress_homo
1247 // )
1248 // );
1249 // }
1250 // for(
1251 // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
1252 // mit!=trans_elastic.setOfBlocks.end();
1253 // mit++
1254 // ) {
1255 // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
1256 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1257 // nitsche_element_trans.getOpPtrVector().push_back(
1258 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1259 // "DISPLACEMENT",elastic_common_data
1260 // )
1261 // );
1262 // nitsche_element_trans.getOpPtrVector().push_back(
1263 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1264 // "MESH_NODE_POSITIONS",elastic_common_data
1265 // )
1266 // );
1267 // nitsche_element_trans.getOpPtrVector().push_back(
1268 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1269 // "POTENTIAL_FIELD",elastic_common_data
1270 // )
1271 // );
1272 // nitsche_element_trans.getOpPtrVector().push_back(
1273 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1274 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,false,false,true
1275 // )
1276 // );
1277 // nitsche_element_trans.getOpPtrVector().push_back(
1278 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral(
1279 // "DISPLACEMENT",nitsche_block_data,nitsche_common_data,
1280 // elastic_block_data,elastic_common_data,stress_homo
1281 // )
1282 // );
1283 // }
1284 // } else {
1285 // for(
1286 // map<int,NonlinearElasticElement::BlockData>::iterator mit = iso_elastic.setOfBlocks.begin();
1287 // mit!=iso_elastic.setOfBlocks.end();
1288 // mit++
1289 // ) {
1290 // NonlinearElasticElement::CommonData &elastic_common_data = iso_elastic.commonData;
1291 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1292 // ave_stress_iso.getOpPtrVector().push_back(
1293 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1294 // "DISPLACEMENT",elastic_common_data
1295 // )
1296 // );
1297 // ave_stress_iso.getOpPtrVector().push_back(
1298 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1299 // "MESH_NODE_POSITIONS",elastic_common_data
1300 // )
1301 // );
1302 // ave_stress_iso.getOpPtrVector().push_back(
1303 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1304 // "DISPLACEMENT",elastic_block_data,elastic_common_data,1,false,false,true
1305 // )
1306 // );
1307 // ave_stress_iso.getOpPtrVector().push_back(
1308 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral(
1309 // "DISPLACEMENT",elastic_block_data,elastic_common_data,stress_homo
1310 // )
1311 // );
1312 // }
1313 // for(
1314 // map<int,NonlinearElasticElement::BlockData>::iterator mit = trans_elastic.setOfBlocks.begin();
1315 // mit!=trans_elastic.setOfBlocks.end();
1316 // mit++
1317 // ) {
1318 // NonlinearElasticElement::CommonData &elastic_common_data = trans_elastic.commonData;
1319 // NonlinearElasticElement::BlockData &elastic_block_data = mit->second;
1320 // ave_stress_trans.getOpPtrVector().push_back(
1321 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1322 // "DISPLACEMENT",elastic_common_data
1323 // )
1324 // );
1325 // ave_stress_trans.getOpPtrVector().push_back(
1326 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1327 // "MESH_NODE_POSITIONS",elastic_common_data
1328 // )
1329 // );
1330 // ave_stress_trans.getOpPtrVector().push_back(
1331 // new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
1332 // "POTENTIAL_FIELD",elastic_common_data
1333 // )
1334 // );
1335 // ave_stress_trans.getOpPtrVector().push_back(
1336 // new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
1337 // "DISPLACEMENT",elastic_block_data,elastic_common_data,2,false,false,true
1338 // )
1339 // );
1340 // ave_stress_trans.getOpPtrVector().push_back(
1341 // new PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral(
1342 // "DISPLACEMENT",elastic_block_data,elastic_common_data,stress_homo
1343 // )
1344 // );
1345 // }
1346 // }
1347 break;
1348 case NOHOMOBC:
1349 SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"*** Boundary conditions not set");
1350 }
1351
1352 struct MyPostProc: public PostProcVolumeOnRefinedMesh {
1353
1354 bool doPreProcess;
1355 bool doPostProcess;
1356
1357 MyPostProc(MoFEM::Interface &m_field):
1358 PostProcVolumeOnRefinedMesh(m_field),
1359 doPreProcess(true),
1360 doPostProcess(true)
1361 {}
1362
1363 void setDoPreProcess() { doPreProcess = true; }
1364 void unSetDoPreProcess() { doPreProcess = false; }
1365 void setDoPostProcess() { doPostProcess = true; }
1366 void unSetDoPostProcess() { doPostProcess = false; }
1367
1368
1369
1370 PetscErrorCode preProcess() {
1371 PetscFunctionBegin;
1372 if(doPreProcess) {
1373 ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
1374 }
1375 PetscFunctionReturn(0);
1376 }
1377 PetscErrorCode postProcess() {
1378 PetscFunctionBegin;
1379 if(doPostProcess) {
1380 ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
1381 }
1382 PetscFunctionReturn(0);
1383 }
1384 };
1385
1386 MyPostProc post_proc(m_field);
1387
1388 ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
1389 ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1390 ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
1391 ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1392 ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
1393 if(trans_iso_blocks) {
1394 ierr = post_proc.addFieldValuesGradientPostProc("POTENTIAL_FIELD"); CHKERRQ(ierr);
1395 }
1396 for(
1397 map<int,NonlinearElasticElement::BlockData>::iterator sit = iso_elastic.setOfBlocks.begin();
1398 sit != iso_elastic.setOfBlocks.end(); sit++
1399 ) {
1400 post_proc.getOpPtrVector().push_back(
1401 new PostPorcStress(
1402 post_proc.postProcMesh,
1403 post_proc.mapGaussPts,
1404 "DISPLACEMENT",
1405 sit->second,
1406 post_proc.commonData,
1407 false
1408 )
1409 );
1410 }
1411 for(
1412 map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
1413 sit != trans_elastic.setOfBlocks.end(); sit++
1414 ) {
1415 post_proc.getOpPtrVector().push_back(
1416 new PostPorcStress(
1417 post_proc.postProcMesh,
1418 post_proc.mapGaussPts,
1419 "DISPLACEMENT",
1420 sit->second,
1421 post_proc.commonData
1422 )
1423 );
1424 }
1425
1426 PetscScalar *avec;
1427 ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
1428 for(int ii = 0;ii<6;ii++) {
1429 ierr = VecZeroEntries(D); CHKERRQ(ierr);
1430 ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
1431 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1432 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1433 ierr = m_field.set_global_ghost_vector(
1434 "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
1435 ); CHKERRQ(ierr);
1436 post_proc.setDoPreProcess();
1437 post_proc.unSetDoPostProcess();
1438 ierr = m_field.loop_finite_elements(
1439 "ELASTIC_MECHANICS","ELASTIC",post_proc
1440 ); CHKERRQ(ierr);
1441 post_proc.unSetDoPreProcess();
1442 post_proc.setDoPostProcess();
1443
1444 if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1445 ierr = m_field.loop_finite_elements(
1446 "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",post_proc
1447 ); CHKERRQ(ierr);
1448 {
1449 ostringstream sss;
1450 sss << "mode_" << homo_bc_names[choise_value] << "_" << ii << ".h5m";
1451 rval = post_proc.postProcMesh.write_file(
1452 sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
1453 ); CHKERRQ_MOAB(rval);
1454 }
1455 ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
1456 if(choise_value == HOMOBCDISP) {
1457 ierr = m_field.loop_finite_elements(
1458 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_disp.getLoopFeRVEBCStress()
1459 ); CHKERRQ(ierr);
1460 }
1461 if(choise_value == HOMOBCTRAC) {
1462 ierr = m_field.loop_finite_elements(
1463 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_trac.getLoopFeRVEBCStress()
1464 ); CHKERRQ(ierr);
1465 }
1466 if(choise_value == HOMOBCPERIODIC) {
1467 ierr = m_field.loop_finite_elements(
1468 "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
1469 ); CHKERRQ(ierr);
1470 }
1471 // if(choise_value == NITSCHE) {
1472 // if(stress_by_boundary_integral) {
1473 // nitsche_element_iso.addToRule = 1;
1474 // nitsche_element_trans.addToRule = 1;
1475 // ierr = m_field.loop_finite_elements(
1476 // "ELASTIC_MECHANICS","ELASTIC",nitsche_element_iso
1477 // ); CHKERRQ(ierr);
1478 // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1479 // ierr = m_field.loop_finite_elements(
1480 // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",nitsche_element_trans
1481 // ); CHKERRQ(ierr);
1482 // } else {
1483 // ave_stress_iso.addToRule = 1;
1484 // ave_stress_trans.addToRule = 1;
1485 // ierr = m_field.loop_finite_elements(
1486 // "ELASTIC_MECHANICS","ELASTIC",ave_stress_iso
1487 // ); CHKERRQ(ierr);
1488 // if(m_field.check_finite_element("TRAN_ISOTROPIC_ELASTIC"))
1489 // ierr = m_field.loop_finite_elements(
1490 // "ELASTIC_MECHANICS","TRAN_ISOTROPIC_ELASTIC",ave_stress_trans
1491 // ); CHKERRQ(ierr);
1492 // }
1493 // }
1494 ierr = PetscOptionsGetReal(
1495 PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
1496 ); CHKERRQ(ierr);
1497 ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
1498 ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
1499 ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
1500 ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1501 ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
1502 for(int jj=0; jj<6; jj++) {
1503 Dmat(jj,ii) = avec[jj];
1504 }
1505 }
1506 ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
1507
1508 PetscPrintf(
1509 PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
1510 );
1511
1512 for(int ii=0; ii<6; ii++) {
1513 PetscPrintf(
1514 PETSC_COMM_WORLD,
1515 "stress %d\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\n",
1516 ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
1517 );
1518 }
1519
1520 //Saving Dmat as a bindary file to use it macro-structure
1521
1522
1523 char output_file_Dmat[255];
1524 ierr = PetscOptionsGetString(PETSC_NULL,"-my_output_file_Dmat",output_file_Dmat,255,&flg); CHKERRQ(ierr);
1525 if(flg) {
1526
1527 //Reading and writing binary files
1528 if(pcomm->rank()==0){
1529 int fd;
1530 PetscViewer view_out;
1531 PetscViewerBinaryOpen(PETSC_COMM_SELF,output_file_Dmat,FILE_MODE_WRITE,&view_out);
1532 PetscViewerBinaryGetDescriptor(view_out,&fd);
1533 PetscBinaryWrite(fd,&Dmat(0,0),36,PETSC_DOUBLE,PETSC_FALSE);
1534 PetscViewerDestroy(&view_out);
1535 }
1536
1537 // MatrixDouble Dmat1;
1538 // Dmat1.resize(6,6); Dmat1.clear();
1539 // if(pcomm->rank()==0){
1540 // int fd;
1541 // PetscViewer view_in;
1542 // PetscViewerBinaryOpen(PETSC_COMM_SELF,output_file_Dmat,FILE_MODE_READ,&view_in);
1543 // PetscViewerBinaryGetDescriptor(view_in,&fd);
1544 // PetscBinaryRead(fd,&Dmat1(0,0),36,PETSC_DOUBLE);
1545 // PetscViewerDestroy(&view_in);
1546 // cout<< "Dmat1 After Reading= "<<Dmat1<<endl;
1547 // }
1548 //
1549 }
1550
1551 //detroy matrices
1552 for(int ii = 0;ii<6;ii++) {
1553 ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
1554 }
1555 ierr = VecDestroy(&D); CHKERRQ(ierr);
1556 ierr = MatDestroy(&Aij); CHKERRQ(ierr);
1557 ierr = KSPDestroy(&solver); CHKERRQ(ierr);
1558 ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
1559
1560 ierr = PetscTime(&v2);CHKERRQ(ierr);
1561 ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
1562
1563 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
1564
1565 }
1566 CATCH_ERRORS;
1567
1568
1569 MoFEM::Core::Finalize();
1570}
BCs_RVELagrange_Disp.hpp
BCs_RVELagrange_Periodic.hpp
BCs_RVELagrange_Trac.hpp
BCs_RVELagrange_Trac_Rigid_Rot.hpp
BCs_RVELagrange_Trac_Rigid_Trans.hpp
Hooke.hpp
Implementation of linear elastic material.
NonLinearElasticElement.hpp
Operators and data structures for non-linear elastic analysis.
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
PostProcStresses.hpp
Post-processing stresses for non-linear analysis.
PostPorcStress
DEPRECATED typedef PostProcStress PostPorcStress
Definition: PostProcStresses.hpp:193
VolumeCalculation.hpp
Operator can be used with any volume element to calculate sum of volumes of all volumes in the set.
main
int main()
Definition: adol-c_atom.cpp:46
a
constexpr double a
Definition: approx_sphere.cpp:30
COL
@ COL
Definition: definitions.h:123
ROW
@ ROW
Definition: definitions.h:123
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
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:541
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:215
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
SIDESET
@ SIDESET
Definition: definitions.h:147
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
MOFEM_IMPOSIBLE_CASE
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
order
int order
Definition: fluid_structure_eigenproblem.cpp:84
eta
constexpr double eta
Definition: free_surface.cpp:96
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::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::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_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::CoreInterface::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
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
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
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
D
double D
Definition: initial_diffusion.cpp:44
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
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
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:152
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
zeta
double zeta
Definition: plastic.cpp:104
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: rve_mechanical.cpp:134
help
static char help[]
Definition: rve_mechanical.cpp:77
tricircumcenter3d_tp
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
HomoBCTypes
HomoBCTypes
Definition: rve_mechanical.cpp:79
NOHOMOBC
@ NOHOMOBC
Definition: rve_mechanical.cpp:84
HOMOBCPERIODIC
@ HOMOBCPERIODIC
Definition: rve_mechanical.cpp:81
NITSCHE
@ NITSCHE
Definition: rve_mechanical.cpp:83
HOMOBCDISP
@ HOMOBCDISP
Definition: rve_mechanical.cpp:80
HOMOBCTRAC
@ HOMOBCTRAC
Definition: rve_mechanical.cpp:82
homo_bc_names
const char * homo_bc_names[]
Definition: rve_mechanical.cpp:87
tetcircumcenter_tp
void tetcircumcenter_tp(double a[3], double b[3], double c[3], double d[3], double circumcenter[3], double *xi, double *eta, double *zeta)
BCs_RVELagrange_Disp
Definition: BCs_RVELagrange_Disp.hpp:18
BCs_RVELagrange_Disp::setOfRVEBC
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
Definition: BCs_RVELagrange_Disp.hpp:56
BCs_RVELagrange_Disp::getLoopFeRVEBCLhs
MyTriFE & getLoopFeRVEBCLhs()
Definition: BCs_RVELagrange_Disp.hpp:35
BCs_RVELagrange_Disp::getLoopFeRVEBCRhs
MyTriFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Disp.hpp:36
BCs_RVELagrange_Disp::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec Stress_Homo)
Definition: BCs_RVELagrange_Disp.hpp:671
BCs_RVELagrange_Disp::getLoopFeRVEBCStress
MyTriFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Disp.hpp:39
BCs_RVELagrange_Disp::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec)
Definition: BCs_RVELagrange_Disp.hpp:580
BCs_RVELagrange_Disp::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
Definition: BCs_RVELagrange_Disp.hpp:58
BCs_RVELagrange_Periodic
Definition: BCs_RVELagrange_Periodic.hpp:19
BCs_RVELagrange_Periodic::getLoopFeRVEBCRhs
MyPrismFE & getLoopFeRVEBCRhs()
Definition: BCs_RVELagrange_Periodic.hpp:39
BCs_RVELagrange_Periodic::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Definition: BCs_RVELagrange_Periodic.hpp:1312
BCs_RVELagrange_Periodic::getLoopFeRVEBCStress
MyPrismFE & getLoopFeRVEBCStress()
Definition: BCs_RVELagrange_Periodic.hpp:42
BCs_RVELagrange_Periodic::addLagrangiangElement
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
Definition: BCs_RVELagrange_Periodic.hpp:52
BCs_RVELagrange_Periodic::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
Definition: BCs_RVELagrange_Periodic.hpp:1176
BCs_RVELagrange_Periodic::getLoopFeRVEBCLhs
MyPrismFE & getLoopFeRVEBCLhs()
Definition: BCs_RVELagrange_Periodic.hpp:38
BCs_RVELagrange_Trac_Rigid_Rot
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:22
BCs_RVELagrange_Trac_Rigid_Rot::setRVEBCsRigidBodyRotOperators
PetscErrorCode setRVEBCsRigidBodyRotOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Definition: BCs_RVELagrange_Trac_Rigid_Rot.hpp:105
BCs_RVELagrange_Trac_Rigid_Trans
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:18
BCs_RVELagrange_Trac_Rigid_Trans::setRVEBCsRigidBodyTranOperators
PetscErrorCode setRVEBCsRigidBodyTranOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
Definition: BCs_RVELagrange_Trac_Rigid_Trans.hpp:109
BCs_RVELagrange_Trac
Definition: BCs_RVELagrange_Trac.hpp:22
BCs_RVELagrange_Trac::setRVEBCsOperators
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
Definition: BCs_RVELagrange_Trac.hpp:383
BCs_RVELagrange_Trac::setRVEBCsHomoStressOperators
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Definition: BCs_RVELagrange_Trac.hpp:448
Composite_xyzcoord
Definition: homogenisation_periodic_atom_test.cpp:72
Face_CenPos_Handle
Definition: homogenisation_periodic_atom_test.cpp:61
Face_CenPos_Handle::zcoord
double zcoord
Definition: homogenisation_periodic_atom_test.cpp:62
Face_CenPos_Handle::Tri_Hand
const EntityHandle Tri_Hand
Definition: homogenisation_periodic_atom_test.cpp:63
Face_CenPos_Handle::Face_CenPos_Handle
Face_CenPos_Handle(double _xcoord, double _ycoord, double _zcoord, const EntityHandle _Tri_Hand)
Definition: rve_mechanical.cpp:102
Face_CenPos_Handle::xcoord
double xcoord
Definition: homogenisation_periodic_atom_test.cpp:62
Face_CenPos_Handle::ycoord
double ycoord
Definition: homogenisation_periodic_atom_test.cpp:62
Hooke
Hook equation.
Definition: Hooke.hpp:21
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_field_meshset
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::DeprecatedCoreInterface::seed_ref_level
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
Definition: DeprecatedCoreInterface.cpp:24
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
Definition: DeprecatedCoreInterface.cpp:18
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
MoFEM::Mat_Elastic_TransIso
Transverse Isotropic material data structure.
Definition: MaterialBlocks.hpp:371
MoFEM::Mat_Elastic_TransIso::data
_data_ data
Definition: MaterialBlocks.hpp:384
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:157
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::SnesMethod::snes_x
Vec & snes_x
state vector
Definition: LoopMethods.hpp:121
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition: LoopMethods.hpp:124
MyPostProc
Definition: plot_base.cpp:34
MyPostProc::postProcess
MoFEMErrorCode postProcess()
Definition: plot_base.cpp:622
MyPostProc::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Definition: plot_base.cpp:397
MyPostProc::preProcess
MoFEMErrorCode preProcess()
Definition: plot_base.cpp:612
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:557
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
NonlinearElasticElement::getLoopFeLhs
MyVolumeFE & getLoopFeLhs()
get lhs volume element
Definition: NonLinearElasticElement.hpp:68
NonlinearElasticElement::feRhs
MyVolumeFE feRhs
calculate right hand side for tetrahedral elements
Definition: NonLinearElasticElement.hpp:65
NonlinearElasticElement::addElement
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
Definition: NonLinearElasticElement.cpp:1120
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
NonlinearElasticElement::setOperators
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
Definition: NonLinearElasticElement.cpp:1153
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::preProcess
MoFEMErrorCode preProcess()
Definition: PostProcOnRefMesh.hpp:695
PostProcTemplateVolumeOnRefinedMesh< MoFEM::VolumeElementForcesAndSourcesCore >::postProcess
MoFEMErrorCode postProcess()
Definition: PostProcOnRefMesh.hpp:796
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:956
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
Generated by Doxygen 1.9.5 and hosted at University of Glasgow