v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
rve_fibre_directions.cpp File Reference
#include <gm_rule.h>
#include <MoFEM.hpp>
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#include <FEMethod_LowLevelStudent.hpp>
#include <FEMethod_UpLevelStudent.hpp>
#include <PostProcVertexMethod.hpp>
#include <PostProcDisplacementAndStrainOnRefindedMesh.hpp>
#include <PotentialFlowFEMethod.hpp>
#include <SurfacePressure.hpp>

Go to the source code of this file.

Functions

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

Variables

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

Function Documentation

◆ main()

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

Populating the Multi-index container with nodes on +ve faces

Getting No. of Fibres and their index to be used for Potential Flow Problem

Adding entities to Field and FE for Potential Flow Problem

Definition at line 48 of file rve_fibre_directions.cpp.

48 {
49
50 try {
51
52 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
53
54 moab::Core mb_instance;
55 moab::Interface& moab = mb_instance;
56 int rank;
57 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
58
59 PetscBool flg = PETSC_TRUE;
60 char mesh_file_name[255];
61 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
62 if(flg != PETSC_TRUE) {
63 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
64 }
65 PetscInt order;
66 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
67 if(flg != PETSC_TRUE) {
68 order = 1;
69 }
70
71 PetscInt mesh_refinement_level;
72 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_mesh_ref_level",&mesh_refinement_level,&flg); CHKERRQ(ierr);
73 if(flg != PETSC_TRUE) {
74 mesh_refinement_level = 0;
75 }
76
77 const char *option;
78 option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
79 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
80 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
81 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
82
83
84 MoFEM::Core core(moab);
85 MoFEM::Interface& m_field = core;
86
87
88 //=======================================================================================================
89 //Seting nodal coordinates on the surface to make sure they are periodic
90 //=======================================================================================================
91
92 Range SurTrisNeg, SurTrisPos;
93 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
94 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
95 ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
96 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
97
98 Range SurNodesNeg,SurNodesPos;
99 rval = moab.get_connectivity(SurTrisNeg,SurNodesNeg,true); CHKERRQ_MOAB(rval);
100 cout<<" All nodes on negative surfaces " << SurNodesNeg.size()<<endl;
101 rval = moab.get_connectivity(SurTrisPos,SurNodesPos,true); CHKERRQ_MOAB(rval);
102 cout<<" All nodes on positive surfaces " << SurNodesPos.size()<<endl;
103
104
105 double roundfact=1e3; double coords_nodes[3];
106 //Populating the Multi-index container with nodes on -ve faces
107 for(Range::iterator nit = SurNodesNeg.begin(); nit!=SurNodesNeg.end(); nit++) {
108 rval = moab.get_coords(&*nit,1,coords_nodes); CHKERRQ_MOAB(rval);
109 //round values to 3 disimal places
110 if(coords_nodes[0]>=0) coords_nodes[0]=double(int(coords_nodes[0]*roundfact+0.5))/roundfact; else coords_nodes[0]=double(int(coords_nodes[0]*roundfact-0.5))/roundfact;
111 if(coords_nodes[1]>=0) coords_nodes[1]=double(int(coords_nodes[1]*roundfact+0.5))/roundfact; else coords_nodes[1]=double(int(coords_nodes[1]*roundfact-0.5))/roundfact;
112 if(coords_nodes[2]>=0) coords_nodes[2]=double(int(coords_nodes[2]*roundfact+0.5))/roundfact; else coords_nodes[2]=double(int(coords_nodes[2]*roundfact-0.5))/roundfact;
113 rval = moab.set_coords(&*nit,1,coords_nodes); CHKERRQ_MOAB(rval);
114 // cout<<" coords_nodes[0]= "<<coords_nodes[0] << " coords_nodes[1]= "<< coords_nodes[1] << " coords_nodes[2]= "<< coords_nodes[2] <<endl;
115 }
116
117 ///Populating the Multi-index container with nodes on +ve faces
118 for(Range::iterator nit = SurNodesPos.begin(); nit!=SurNodesPos.end(); nit++) {
119 rval = moab.get_coords(&*nit,1,coords_nodes); CHKERRQ_MOAB(rval);
120 //round values to 3 disimal places
121 if(coords_nodes[0]>=0) coords_nodes[0]=double(int(coords_nodes[0]*roundfact+0.5))/roundfact; else coords_nodes[0]=double(int(coords_nodes[0]*roundfact-0.5))/roundfact;
122 if(coords_nodes[1]>=0) coords_nodes[1]=double(int(coords_nodes[1]*roundfact+0.5))/roundfact; else coords_nodes[1]=double(int(coords_nodes[1]*roundfact-0.5))/roundfact;
123 if(coords_nodes[2]>=0) coords_nodes[2]=double(int(coords_nodes[2]*roundfact+0.5))/roundfact; else coords_nodes[2]=double(int(coords_nodes[2]*roundfact-0.5))/roundfact;
124 rval = moab.set_coords(&*nit,1,coords_nodes); CHKERRQ_MOAB(rval);
125 // cout<<" coords_nodes[0]= "<<coords_nodes[0] << " coords_nodes[1]= "<< coords_nodes[1] << " coords_nodes[2]= "<< coords_nodes[2] <<endl;
126 }
127 //=======================================================================================================
128
129 //add fields
130 ierr = m_field.add_field("POTENTIAL_FIELD",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
131 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
132
133 ///Getting No. of Fibres and their index to be used for Potential Flow Problem
134 int noOfFibres=0;
135 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
136
137 std::size_t found=it->getName().find("PotentialFlow");
138 if (found==std::string::npos) continue;
139 noOfFibres += 1;
140 }
141 cout<<"No. of Fibres for Potential Flow : "<<noOfFibres<<endl;
142
143 vector<int> fibreList(noOfFibres,0);
144 int aaa=0;
145
146 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
147
148 std::size_t interfaceFound=it->getName().find("PotentialFlow_");
149 if (interfaceFound==std::string::npos) continue;
150
151 std::string str2 = it->getName().substr (14,50);
152
153 fibreList[aaa] = atoi(str2.c_str());
154 aaa += 1;
155 }
156
157 //************************************************************//
158
159 for (int cc = 0; cc < noOfFibres; cc++) {
160
161 ostringstream sss, rrr;
162
163 //add finite elements
164 sss << "POTENTIAL_ELEM" << fibreList[cc];
165 cout<<sss.str().c_str()<<endl;
166 ierr = m_field.add_finite_element( sss.str().c_str() ); CHKERRQ(ierr);
167 ierr = m_field.modify_finite_element_add_field_row( sss.str().c_str() ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
168 ierr = m_field.modify_finite_element_add_field_col( sss.str().c_str() ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
169 ierr = m_field.modify_finite_element_add_field_data( sss.str().c_str() ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
170 ierr = m_field.modify_finite_element_add_field_data( sss.str().c_str() ,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
171
172 //add problems
173 rrr << "POTENTIAL_PROBLEM" << fibreList[cc];
174 ierr = m_field.add_problem( rrr.str().c_str() ); CHKERRQ(ierr);
175 //define problems and finite elements
176 ierr = m_field.modify_problem_add_finite_element( rrr.str().c_str() , sss.str().c_str() ); CHKERRQ(ierr);
177
178 }
179
180 Tag th_meshset_info;
181 int def_meshset_info[2] = {0,0};
182 rval = moab.tag_get_handle(
183 "MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info
184 );
185
186 int meshset_data[2];
187 EntityHandle root = moab.get_root_set();
188 rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
189
190 ierr = m_field.seed_ref_level_3D(0,BitRefLevel().set(0)); CHKERRQ(ierr);
191 vector<BitRefLevel> bit_levels;
192 bit_levels.push_back(BitRefLevel().set(0));
193
194 //MESH REFINEMENT
195 int ll = 1;
196 //End of refinement, save level of refinement
197 int meshset_data_root[2]={ll,0};
198 rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data_root); CHKERRQ_MOAB(rval);
199
200 /******************TETS TO MESHSET AND SAVING TETS ENTITIES******************/
201 EntityHandle out_tet_meshset;
202 rval = moab.create_meshset(MESHSET_SET,out_tet_meshset); CHKERRQ_MOAB(rval);
203 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBTET,out_tet_meshset); CHKERRQ(ierr);
204 rval = moab.write_file("out_tets.vtk","VTK","",&out_tet_meshset,1); CHKERRQ_MOAB(rval);
205 /*******************************************************/
206
207 EntityHandle out_meshset;
208 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
209 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
210 rval = moab.write_file("out_all_mesh.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
211 Range LatestRefinedTets;
212 rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
213
214
215 BitRefLevel problem_bit_level = bit_levels.back();
216
217
218 ///Adding entities to Field and FE for Potential Flow Problem
219 for (int cc = 0; cc < noOfFibres; cc++) {
220
221 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
222
223 // std::size_t found=it->getName().find("PotentialFlow");
224 // if (found==std::string::npos) continue;
225 // cout<<it->getName()<<endl;
226
227 ostringstream sss,rrr;
228 //set problem level
229 sss << "POTENTIAL_ELEM" << fibreList[cc];
230 rrr << "PotentialFlow_" << fibreList[cc];
231
232 if(it->getName() == rrr.str().c_str() ) {
233 Range TetsInBlock;
234 rval = moab.get_entities_by_type(it->meshset, MBTET,TetsInBlock,true); CHKERRQ_MOAB(rval);
235 Range block_rope_bit_level = intersect(LatestRefinedTets,TetsInBlock);
236
237 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"POTENTIAL_FIELD"); CHKERRQ(ierr);
238 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
239
240 //add finite elements entities
241 ierr = m_field.add_ents_to_finite_element_by_type(block_rope_bit_level,MBTET, sss.str().c_str()); CHKERRQ(ierr);
242
243 }
244 }
245 }
246
247 ierr = m_field.set_field_order(0,MBVERTEX,"POTENTIAL_FIELD",1); CHKERRQ(ierr);
248 ierr = m_field.set_field_order(0,MBEDGE,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
249 ierr = m_field.set_field_order(0,MBTRI,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
250 ierr = m_field.set_field_order(0,MBTET,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
251
252 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
253 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
254 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
255 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
256
257
258 struct myMetaNeummanForces{
259
260 static PetscErrorCode addNeumannFluxBCElements(
261 MoFEM::Interface &m_field,
262 const string problem_name,
263 const string field_name,
264 const int fibre_id,
265 const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
266
267 PetscFunctionBegin;
268
269
270
271 ostringstream sss,rrr,ppp;
272 ppp << "PressureIO_" << fibre_id <<"_1";
273 sss << "PressureIO_" << fibre_id <<"_2";
274 rrr << "FLUX_FE" <<fibre_id;
275 ierr = m_field.add_finite_element( rrr.str().c_str() ,MF_ZERO); CHKERRQ(ierr);
276 ierr = m_field.modify_finite_element_add_field_row(rrr.str().c_str(),field_name); CHKERRQ(ierr);
277 ierr = m_field.modify_finite_element_add_field_col(rrr.str().c_str(),field_name); CHKERRQ(ierr);
278 ierr = m_field.modify_finite_element_add_field_data(rrr.str().c_str(),field_name); CHKERRQ(ierr);
279 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,SIDESET|PRESSURESET,it)) {
280
281 // std::size_t PressureFound=it->getName().find(sss.str().c_str());
282 // if (PressureFound==std::string::npos) continue;
283 if (ppp.str().c_str()==it->getName() || sss.str().c_str()==it->getName()){
284
285 if(m_field.check_field(mesh_nodals_positions)) {
286 ierr = m_field.modify_finite_element_add_field_data( rrr.str().c_str() ,mesh_nodals_positions); CHKERRQ(ierr);
287 }
288 ierr = m_field.modify_problem_add_finite_element(problem_name, rrr.str().c_str() ); CHKERRQ(ierr);
289 Range tris;
290 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
291 ierr = m_field.add_ents_to_finite_element_by_type(tris,MBTRI, rrr.str().c_str() ); CHKERRQ(ierr);
292 }
293 }
294
295 PetscFunctionReturn(0);
296 }
297
298 static PetscErrorCode setNeumannFluxFiniteElementOperators(
299 MoFEM::Interface &m_field,
300 boost::ptr_map<string,NeummanForcesSurface> &neumann_forces,
301 Vec &F,const string field_name,const int fibre_id, const string mesh_nodals_positions = "MESH_NODE_POSITIONS") {
302
303 PetscFunctionBegin;
304
305 string fe_name;
306 // fe_name = "FLUX_FE";
307 ostringstream sss,rrr,ppp;
308 ppp << "PressureIO_" << fibre_id<<"_1";
309 sss << "PressureIO_" << fibre_id<<"_2";
310 rrr << "FLUX_FE" <<fibre_id;
311 fe_name = rrr.str().c_str();
312 neumann_forces.insert(fe_name,new NeummanForcesSurface(m_field));
313 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,SIDESET|PRESSURESET,it)) {
314
315 // std::size_t PressureFound=it->getName().find(sss.str().c_str());
316 // if (PressureFound==std::string::npos) continue;
317 if (ppp.str().c_str()==it->getName() || sss.str().c_str()==it->getName()){
318 bool ho_geometry = m_field.check_field(mesh_nodals_positions);
319 ierr = neumann_forces.at(fe_name).addFlux(field_name,F,it->getMeshsetId(),ho_geometry); CHKERRQ(ierr);
320 /*pressure_cubit_bc_data data;
321 ierr = it->getBcDataStructure(data); CHKERRQ(ierr);
322 my_split << *it << endl;
323 my_split << data << endl;*/
324 }
325 }
326 PetscFunctionReturn(0);
327 }
328
329 };
330
331 for (int cc = 0; cc < noOfFibres; cc++) {
332 ostringstream sss;
333 sss << "POTENTIAL_PROBLEM" << fibreList[cc];
334
335 //flux boundary conditions
336 ierr = myMetaNeummanForces::addNeumannFluxBCElements(m_field,sss.str().c_str(),"POTENTIAL_FIELD",fibreList[cc]); CHKERRQ(ierr);
337 //set problem level
338 ierr = m_field.modify_problem_ref_level_add_bit( sss.str().c_str() ,problem_bit_level); CHKERRQ(ierr);
339 }
340
341 //build fields
342 ierr = m_field.build_fields(); CHKERRQ(ierr);
343 //build finite elements
344 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
345 //build adjacencies
346 ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
347 //build problem
348 ierr = m_field.build_problems(); CHKERRQ(ierr);
349
350
351 //partition problems
352 for (int cc = 0; cc < noOfFibres; cc++) {
353 ostringstream sss;
354 sss << "POTENTIAL_PROBLEM" << fibreList[cc];
355 ierr = m_field.partition_problem( sss.str().c_str() ); CHKERRQ(ierr);
356 ierr = m_field.partition_finite_elements( sss.str().c_str() ); CHKERRQ(ierr);
357 ierr = m_field.partition_ghost_dofs( sss.str().c_str() ); CHKERRQ(ierr);
358 }
359
360 //print bcs
361 ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
362 ierr = m_field.print_cubit_pressure_set(); CHKERRQ(ierr);
363
364 //create matrices and vectors
365 vector<Vec> F(noOfFibres);
366 vector<Vec> D(noOfFibres);
367 vector<Mat> A(noOfFibres);
368
369 for (int cc = 0; cc < noOfFibres; cc++) {
370 ostringstream sss,rrr,ttt;
371 sss << "POTENTIAL_PROBLEM" << fibreList[cc];
372 rrr << "POTENTIAL_ELEM" << fibreList[cc];
373 ttt << "ZeroPressure_" << fibreList[cc];
374 ierr = m_field.VecCreateGhost( sss.str().c_str() ,ROW,&F[cc]); CHKERRQ(ierr);
375 ierr = m_field.VecCreateGhost( sss.str().c_str() ,COL,&D[cc]); CHKERRQ(ierr);
376 ierr = m_field.MatCreateMPIAIJWithArrays( sss.str().c_str() ,&A[cc]); CHKERRQ(ierr);
377
378 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
379 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
380
381 //get nodes and other entities to fix
382 Range fix_nodes;
383 for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,NODESET|UNKNOWNNAME,it)) {
384 // std::size_t zeroPressureFound=it->getName().find(ttt.str().c_str());
385 // if (zeroPressureFound==std::string::npos) continue;
386 if (ttt.str().c_str()==it->getName()){
387 rval = moab.get_entities_by_type(it->meshset,MBVERTEX,fix_nodes,true); CHKERRQ_MOAB(rval);
388 Range edges;
389 rval = moab.get_entities_by_type(it->meshset,MBEDGE,edges,true); CHKERRQ_MOAB(rval);
390 Range tris;
391 rval = moab.get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
392 Range adj;
393 rval = moab.get_connectivity(tris,adj,true); CHKERRQ_MOAB(rval);
394 fix_nodes.insert(adj.begin(),adj.end());
395 rval = moab.get_connectivity(edges,adj,true); CHKERRQ_MOAB(rval);
396 fix_nodes.insert(adj.begin(),adj.end());
397 rval = moab.get_adjacencies(tris,1,false,edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
398 }
399 }
400 FixBcAtEntities fix_dofs(m_field,"POTENTIAL_FIELD",A[cc],D[cc],F[cc],fix_nodes);
401 //initialize data structure
402 ierr = m_field.problem_basic_method_preProcess( sss.str().c_str() ,fix_dofs); CHKERRQ(ierr);
403
404 //neuman flux bc elements
405 boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
406 ierr = myMetaNeummanForces::setNeumannFluxFiniteElementOperators(m_field,neumann_forces,F[cc],"POTENTIAL_FIELD",fibreList[cc]); CHKERRQ(ierr);
407 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
408 for(;mit!=neumann_forces.end();mit++) {
409 ierr = m_field.loop_finite_elements( sss.str().c_str() ,mit->first,mit->second->getLoopFe()); CHKERRQ(ierr);
410 }
411
412 LaplacianElem elem(m_field,A[cc],F[cc]);
413
414 ierr = MatZeroEntries(A[cc]); CHKERRQ(ierr);
415 ierr = m_field.loop_finite_elements( sss.str().c_str() , rrr.str().c_str() ,elem); CHKERRQ(ierr);
416
417 //post proces fix boundary conditiond
418 ierr = m_field.problem_basic_method_postProcess( sss.str().c_str() ,fix_dofs); CHKERRQ(ierr);
419
420 ierr = MatAssemblyBegin(A[cc],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
421 ierr = MatAssemblyEnd(A[cc],MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
422
423 ierr = VecAssemblyBegin(F[cc]); CHKERRQ(ierr);
424 ierr = VecAssemblyEnd(F[cc]); CHKERRQ(ierr);
425
426 // VecView(F[0],PETSC_VIEWER_STDOUT_WORLD);
427 //set matrix possitives define and symetric for cholesky and icc preceonditionser
428 ierr = MatSetOption(A[cc],MAT_SPD,PETSC_TRUE); CHKERRQ(ierr);
429
430 KSP solver;
431 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
432 ierr = KSPSetOperators(solver,A[cc],A[cc]); CHKERRQ(ierr);
433 ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
434 ierr = KSPSetUp(solver); CHKERRQ(ierr);
435
436 ierr = KSPSolve(solver,F[cc],D[cc]); CHKERRQ(ierr);
437 ierr = VecGhostUpdateBegin(D[cc],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
438 ierr = VecGhostUpdateEnd(D[cc],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
439 ierr = m_field.set_global_ghost_vector( sss.str().c_str() ,ROW,D[cc],INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
440
441 ierr = KSPDestroy(&solver); CHKERRQ(ierr);
442 ierr = VecDestroy(&F[cc]); CHKERRQ(ierr);
443 ierr = VecDestroy(&D[cc]); CHKERRQ(ierr);
444 ierr = MatDestroy(&A[cc]); CHKERRQ(ierr);
445 }
446
447 Tag th_phi;
448 double def_val = 0;
449 rval = moab.tag_get_handle("PHI",1,MB_TYPE_DOUBLE,th_phi,MB_TAG_CREAT|MB_TAG_SPARSE,&def_val); CHKERRQ_MOAB(rval);
450 for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(m_field,"POTENTIAL_FIELD",dof)) {
451 EntityHandle ent = dof->get()->getEnt();
452 double val = dof->get()->getFieldData();
453 rval = moab.tag_set_data(th_phi,&ent,1,&val); CHKERRQ_MOAB(rval);
454 }
455
456 ProjectionFieldOn10NodeTet ent_method_phi_on_10nodeTet(m_field,"POTENTIAL_FIELD",true,false,"PHI");
457 ierr = m_field.loop_dofs("POTENTIAL_FIELD",ent_method_phi_on_10nodeTet); CHKERRQ(ierr);
458 ent_method_phi_on_10nodeTet.setNodes = false;
459 ierr = m_field.loop_dofs("POTENTIAL_FIELD",ent_method_phi_on_10nodeTet); CHKERRQ(ierr);
460
461 if(pcomm->rank()==0) {
462 rval = moab.write_file("solution_RVE.h5m"); CHKERRQ_MOAB(rval);
463 }
464
465 EntityHandle out_meshset1;
466 rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
467 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels[0],BitRefLevel().set(),MBTET,out_meshset1); CHKERRQ(ierr);
468 rval = moab.write_file("solution2.vtk","VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
469
470 if(pcomm->rank()==0) {
471 EntityHandle out_meshset;
472 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
473
474 for (int cc = 0; cc < noOfFibres; cc++) {
475 EntityHandle out_meshset1;
476 rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
477 ostringstream sss,rrr,ttt;
478 sss << "POTENTIAL_PROBLEM" << fibreList[cc];
479 rrr << "POTENTIAL_ELEM" << fibreList[cc];
480 ttt << "out_potential_flow" << fibreList[cc] <<".vtk";
481 ierr = m_field.get_problem_finite_elements_entities( sss.str().c_str() , rrr.str().c_str() ,out_meshset); CHKERRQ(ierr);
482 ierr = m_field.get_problem_finite_elements_entities( sss.str().c_str() , rrr.str().c_str() ,out_meshset1); CHKERRQ(ierr);
483
484 rval = moab.write_file( ttt.str().c_str() ,"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
485 rval = moab.delete_entities(&out_meshset1,1); CHKERRQ_MOAB(rval);
486 }
487
488 rval = moab.write_file("out_potential_flow.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
489 rval = moab.delete_entities(&out_meshset,1); CHKERRQ_MOAB(rval);
490 }
491
492 MoFEM::Core::Finalize();
493
494 }
495 CATCH_ERRORS;
496
497}
FixBcAtEntities
DEPRECATED typedef DirichletFixFieldAtEntitiesBc FixBcAtEntities
Definition DirichletBC.hpp:331
NeummanForcesSurface
DEPRECATED typedef NeumannForcesSurface NeummanForcesSurface
Definition SurfacePressure.hpp:1321
COL
@ COL
Definition definitions.h:136
ROW
@ ROW
Definition definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MF_ZERO
@ MF_ZERO
Definition definitions.h:111
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
H1
@ H1
continuous field
Definition definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
PRESSURESET
@ PRESSURESET
Definition definitions.h:165
NODESET
@ NODESET
Definition definitions.h:159
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
F
@ F
Definition free_surface.cpp:394
_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_
#define _IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
Definition Interface.hpp:1878
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::problem_basic_method_postProcess
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
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
MoFEM::CoreInterface::get_problem_finite_elements_entities
virtual MoFEMErrorCode get_problem_finite_elements_entities(const std::string name, const std::string &fe_name, const EntityHandle meshset)=0
add finite elements to the meshset
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEM::CoreInterface::modify_problem_add_finite_element
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
D
double D
Definition initial_diffusion.cpp:44
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
A
constexpr AssemblyType A
Definition operators_tests.cpp:30
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
help
static char help[]
Definition rve_fibre_directions.cpp:46
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::problem_basic_method_preProcess
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
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::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
Definition DeprecatedCoreInterface.cpp:18
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24

Variable Documentation

◆ help

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

Definition at line 46 of file rve_fibre_directions.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow