v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
rve_fibre_directions_interface.cpp File Reference
#include <gm_rule.h>
#include <MoFEM.hpp>
#include <boost/math/special_functions/round.hpp>
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#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[] )

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

Solve problem potential flow problem for fibres one by one

Definition at line 46 of file rve_fibre_directions_interface.cpp.

46 {
47
48 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
49
50 try {
51
52 moab::Core mb_instance;
53 moab::Interface& moab = mb_instance;
54 int rank;
55 MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
56
57 PetscBool flg = PETSC_TRUE;
58 char mesh_file_name[255];
59 ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
60 if(flg != PETSC_TRUE) {
61 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
62 }
63 PetscInt order;
64 ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
65 if(flg != PETSC_TRUE) {
66 order = 1;
67 }
68// cout<<" order "<<order<<endl;
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
78 const char *option;
79 option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
80 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
81 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
82 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
83
84
85 MoFEM::Core core(moab);
86 MoFEM::Interface& m_field = core;
87 PrismInterface *interface_ptr;
88 ierr = m_field.query_interface(interface_ptr); CHKERRQ(ierr);
89
90 //=======================================================================================================
91 //Seting nodal coordinates on the surface to make sure they are periodic
92 //=======================================================================================================
93
94 Range SurTrisNeg, SurTrisPos;
95 ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
96 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
97 ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
98 ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
99
100 Range SurNodesNeg,SurNodesPos;
101 rval = moab.get_connectivity(SurTrisNeg,SurNodesNeg,true); CHKERRQ_MOAB(rval);
102 cout<<" All nodes on negative surfaces " << SurNodesNeg.size()<<endl;
103 rval = moab.get_connectivity(SurTrisPos,SurNodesPos,true); CHKERRQ_MOAB(rval);
104 cout<<" All nodes on positive surfaces " << SurNodesPos.size()<<endl;
105
106
107 double roundfact=1e3; double coords_nodes[3];
108 for(Range::iterator nit = SurNodesNeg.begin(); nit!=SurNodesNeg.end(); nit++) {
109 rval = moab.get_coords(&*nit,1,coords_nodes); CHKERRQ_MOAB(rval);
110 //round values to 3 disimal places
111 coords_nodes[0]=round(coords_nodes[0]*roundfact)/roundfact;
112 coords_nodes[1]=round(coords_nodes[1]*roundfact)/roundfact;
113 coords_nodes[2]=round(coords_nodes[2]*roundfact)/roundfact;
114 rval = moab.set_coords(&*nit,1,coords_nodes); CHKERRQ_MOAB(rval);
115// cout<<" coords_nodes[0]= "<<coords_nodes[0] << " coords_nodes[1]= "<< coords_nodes[1] << " coords_nodes[2]= "<< coords_nodes[2] <<endl;
116 }
117
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
122 coords_nodes[0]=round(coords_nodes[0]*roundfact)/roundfact;
123 coords_nodes[1]=round(coords_nodes[1]*roundfact)/roundfact;
124 coords_nodes[2]=round(coords_nodes[2]*roundfact)/roundfact;
125 rval = moab.set_coords(&*nit,1,coords_nodes); CHKERRQ_MOAB(rval);
126// cout<<" coords_nodes[0]= "<<coords_nodes[0] << " coords_nodes[1]= "<< coords_nodes[1] << " coords_nodes[2]= "<< coords_nodes[2] <<endl;
127 }
128 //=======================================================================================================
129
130// string wait;
131// cin>>wait;
132
133 //add fields
134 ierr = m_field.add_field("POTENTIAL_FIELD",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
135 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
136
137 ///Getting No. of Fibres and their index to be used for Potential Flow Problem
138 int noOfFibres=0;
140
141 std::size_t found=it->getName().find("PotentialFlow");
142 if (found==std::string::npos) continue;
143 noOfFibres += 1;
144 }
145 cout<<"No. of Fibres for Potential Flow : "<<noOfFibres<<endl;
146
147 vector<int> fibreList(noOfFibres,0);
148 int aaa=0;
149
151
152 std::size_t interfaceFound=it->getName().find("PotentialFlow_");
153 if (interfaceFound==std::string::npos) continue;
154
155 std::string str2 = it->getName().substr (14,50);
156
157 fibreList[aaa] = atoi(str2.c_str()); //atoi converet string to integer
158// cout<<"atoi(str2.c_str()) = " <<atoi(str2.c_str())<<endl;
159 aaa += 1;
160 }
161
162 //add finite elements
163 ierr = m_field.add_finite_element( "POTENTIAL_ELEM" ); CHKERRQ(ierr);
164 ierr = m_field.modify_finite_element_add_field_row("POTENTIAL_ELEM" ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
165 ierr = m_field.modify_finite_element_add_field_col("POTENTIAL_ELEM" ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
166 ierr = m_field.modify_finite_element_add_field_data("POTENTIAL_ELEM" ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
167 ierr = m_field.modify_finite_element_add_field_data("POTENTIAL_ELEM" ,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
168
169 //add problems
170 ierr = m_field.add_problem( "POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
171 //define problems and finite elements
172 ierr = m_field.modify_problem_add_finite_element( "POTENTIAL_PROBLEM" , "POTENTIAL_ELEM" ); CHKERRQ(ierr);
173
174 Tag th_meshset_info;
175 int def_meshset_info[2] = {0,0};
176 rval = moab.tag_get_handle("MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info);
177
178 int meshset_data[2];
179 EntityHandle root = moab.get_root_set();
180 rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
181
182 ierr = m_field.seed_ref_level_3D(0,BitRefLevel().set(0)); CHKERRQ(ierr);
183 vector<BitRefLevel> bit_levels;
184 bit_levels.push_back(BitRefLevel().set(0));
185
186 //MESH REFINEMENT
187 int ll = 1;
188
189
190 //*****INTERFACE INSERTION******
192
193 std::size_t interfaceFound=cit->getName().find("interface");
194 if (interfaceFound==std::string::npos) continue;
195
196 ierr = PetscPrintf(PETSC_COMM_WORLD,"Insert Interface %d\n",cit->getMeshsetId()); CHKERRQ(ierr);
197 EntityHandle cubit_meshset = cit->getMeshset();
198 int meshset_data_int[2];
199 rval = moab.tag_get_data(th_meshset_info,&cubit_meshset,1,meshset_data_int); CHKERRQ_MOAB(rval);
200 if(meshset_data_int[1]==0){
201 //get tet enties form back bit_level
202 EntityHandle ref_level_meshset = 0;
203 rval = moab.create_meshset(MESHSET_SET,ref_level_meshset); CHKERRQ_MOAB(rval);
204 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBTET,ref_level_meshset); CHKERRQ(ierr);
205 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBPRISM,ref_level_meshset); CHKERRQ(ierr);
206 Range ref_level_tets;
207 rval = moab.get_entities_by_handle(ref_level_meshset,ref_level_tets,true); CHKERRQ_MOAB(rval);
208 //get faces and test to split
209 ierr = interface_ptr->getSides(cubit_meshset,bit_levels.back(),true); CHKERRQ(ierr);
210 //set new bit level
211 bit_levels.push_back(BitRefLevel().set(ll++));
212 //split faces and edges
213 ierr = interface_ptr->splitSides(ref_level_meshset,bit_levels.back(),cubit_meshset,true,true,0); CHKERRQ(ierr);
214 //clean meshsets
215 rval = moab.delete_entities(&ref_level_meshset,1); CHKERRQ_MOAB(rval);
216 int meshset_data_ins[2] = {ll,1};
217 rval = moab.tag_set_data(th_meshset_info,&cubit_meshset,1,meshset_data_ins); CHKERRQ_MOAB(rval);
218 }
219 //update cubit meshsets
220 for(_IT_CUBITMESHSETS_FOR_LOOP_(m_field,ciit)) {
221 EntityHandle cubit_meshset = ciit->meshset;
222 ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBVERTEX,true); CHKERRQ(ierr);
223 ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBEDGE,true); CHKERRQ(ierr);
224 ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBTRI,true); CHKERRQ(ierr);
225 ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBTET,true); CHKERRQ(ierr);
226 }
227 }
228
229
230 //End of refinement, save level of refinement
231 int meshset_data_root[2]={ll,0};
232 rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data_root); CHKERRQ_MOAB(rval);
233
234 /******************TETS TO MESHSET AND SAVING TETS ENTITIES******************/
235 EntityHandle out_tet_meshset;
236 rval = moab.create_meshset(MESHSET_SET,out_tet_meshset); CHKERRQ_MOAB(rval);
237 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBTET,out_tet_meshset); CHKERRQ(ierr);
238 rval = moab.write_file("out_tets.vtk","VTK","",&out_tet_meshset,1); CHKERRQ_MOAB(rval);
239 /*******************************************************/
240
241 /******************PRISMS TO MESHSET AND SAVING PRISMS ENTITIES******************/
242 EntityHandle out_meshset_prism;
243 rval = moab.create_meshset(MESHSET_SET,out_meshset_prism); CHKERRQ_MOAB(rval);
244 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBPRISM,out_meshset_prism); CHKERRQ(ierr);
245 rval = moab.write_file("out_prism.vtk","VTK","",&out_meshset_prism,1); CHKERRQ_MOAB(rval);
246 /*******************************************************/
247
248 EntityHandle out_meshset;
249 rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
250 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
251 rval = moab.write_file("out_all_mesh.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
252
253 Range LatestRefinedTets;
254 rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
255
256 Range LatestRefinedTris;
257 rval = moab.get_entities_by_type(out_meshset, MBTRI,LatestRefinedTris,true); CHKERRQ_MOAB(rval);
258
259
260 BitRefLevel problem_bit_level = bit_levels.back();
261
262 //add entities to field
263 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"POTENTIAL_FIELD"); CHKERRQ(ierr);
264 ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
265
266
267 ierr = m_field.set_field_order(0,MBVERTEX,"POTENTIAL_FIELD",1); CHKERRQ(ierr);
268 ierr = m_field.set_field_order(0,MBEDGE,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
269 ierr = m_field.set_field_order(0,MBTRI,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
270 ierr = m_field.set_field_order(0,MBTET,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
271
272 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
273 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
274 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
275 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
276
277 //define flux element
278 ierr = m_field.add_finite_element("FLUX_FE" ,MF_ZERO); CHKERRQ(ierr);
279 ierr = m_field.modify_finite_element_add_field_row("FLUX_FE","POTENTIAL_FIELD"); CHKERRQ(ierr);
280 ierr = m_field.modify_finite_element_add_field_col("FLUX_FE","POTENTIAL_FIELD"); CHKERRQ(ierr);
281 ierr = m_field.modify_finite_element_add_field_data("FLUX_FE","POTENTIAL_FIELD"); CHKERRQ(ierr);
282 ierr = m_field.modify_finite_element_add_field_data( "FLUX_FE" ,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
283 ierr = m_field.modify_problem_add_finite_element("POTENTIAL_PROBLEM", "FLUX_FE" ); CHKERRQ(ierr);
284
285
286 //create matrices and vectors
287 Vec F;
288 Vec D;
289 Mat A;
290
291 EntityHandle out_meshset_fibres;
292 if(pcomm->rank()==0) {
293 rval = moab.create_meshset(MESHSET_SET,out_meshset_fibres); CHKERRQ_MOAB(rval);
294 }
295
296 ///Solve problem potential flow problem for fibres one by one
297 for (int cc = 0; cc < noOfFibres; cc++) {
298 ostringstream rrr, ppp1, ppp2;
299 rrr << "PotentialFlow_" << fibreList[cc];
300 Range new_tets;
302 if(it->getName() == rrr.str().c_str() ) {
303 //set problem level
304 Range TetsInBlock;
305 rval = moab.get_entities_by_type(it->meshset, MBTET,TetsInBlock,true); CHKERRQ_MOAB(rval);
306 new_tets = intersect(LatestRefinedTets,TetsInBlock);
307 //add finite elements entities
308 ierr = m_field.add_ents_to_finite_element_by_type(new_tets,MBTET,"POTENTIAL_ELEM"); CHKERRQ(ierr);
309 }
310 }
311
312 Range new_tris_1, new_tris_2;
313 //flux boundary conditions
314 ppp1 << "PressureIO_" << fibreList[cc] <<"_1";
315 ppp2 << "PressureIO_" << fibreList[cc] <<"_2";
317 if (ppp1.str().c_str()==it->getName()){
318 Range tris;
319 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
320 new_tris_1 = intersect(LatestRefinedTris,tris);
321 ierr = m_field.add_ents_to_finite_element_by_type(new_tris_1,MBTRI, "FLUX_FE" ); CHKERRQ(ierr);
322 }
323
324
325 if (ppp2.str().c_str()==it->getName()){
326 Range tris;
327 rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
328 new_tris_2 = intersect(LatestRefinedTris,tris);
329 ierr = m_field.add_ents_to_finite_element_by_type(new_tris_2,MBTRI, "FLUX_FE" ); CHKERRQ(ierr);
330 }
331 }
332
333 //set problem level
334 ierr = m_field.modify_problem_ref_level_add_bit("POTENTIAL_PROBLEM" ,problem_bit_level); CHKERRQ(ierr);
335
336 //build fields
337 ierr = m_field.build_fields(); CHKERRQ(ierr);
338 //build finite elements
339 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
340 //build adjacencies
341 ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
342 //build problem
343 ierr = m_field.build_problems(); CHKERRQ(ierr);
344
345 ierr = m_field.partition_problem("POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
346 ierr = m_field.partition_finite_elements( "POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
347 ierr = m_field.partition_ghost_dofs( "POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
348
349 ierr = m_field.VecCreateGhost("POTENTIAL_PROBLEM" ,ROW,&F); CHKERRQ(ierr);
350 ierr = m_field.VecCreateGhost("POTENTIAL_PROBLEM" ,COL,&D); CHKERRQ(ierr);
351 ierr = m_field.MatCreateMPIAIJWithArrays("POTENTIAL_PROBLEM" ,&A); CHKERRQ(ierr);
352
353 Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
354 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
355
356 ostringstream zero_pressure;
357 zero_pressure << "ZeroPressure_" << fibreList[cc];
358 //get nodes and other entities to fix
359 Range fix_nodes;
361 if (zero_pressure.str().c_str()==it->getName()){
362 rval = moab.get_entities_by_type(it->meshset,MBVERTEX,fix_nodes,true); CHKERRQ_MOAB(rval);
363 Range edges;
364 rval = moab.get_entities_by_type(it->meshset,MBEDGE,edges,true); CHKERRQ_MOAB(rval);
365 Range tris;
366 rval = moab.get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
367 Range adj;
368 rval = moab.get_connectivity(tris,adj,true); CHKERRQ_MOAB(rval);
369 fix_nodes.insert(adj.begin(),adj.end());
370 rval = moab.get_connectivity(edges,adj,true); CHKERRQ_MOAB(rval);
371 fix_nodes.insert(adj.begin(),adj.end());
372 rval = moab.get_adjacencies(tris,1,false,edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
373 }
374 }
375 FixBcAtEntities fix_dofs(m_field,"POTENTIAL_FIELD",A,D,F,fix_nodes);
376 //initialize data structure
377 ierr = m_field.problem_basic_method_preProcess("POTENTIAL_PROBLEM",fix_dofs); CHKERRQ(ierr);
378
379 ierr = VecZeroEntries(F); CHKERRQ(ierr);
380 ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
381 ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
382 ierr = MatZeroEntries(A); CHKERRQ(ierr);
383
384 //neuman flux bc elements
385 //surface forces
386 boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
387 string fe_name_str = "FLUX_FE";
388 neumann_forces.insert(fe_name_str,new NeummanForcesSurface(m_field));
390 if (ppp1.str().c_str()==it->getName() || ppp2.str().c_str()==it->getName()){
391 bool ho_geometry = m_field.check_field("MESH_NODE_POSITIONS");
392 ierr = neumann_forces.at("FLUX_FE").addFlux("POTENTIAL_FIELD",F,it->getMeshsetId(),ho_geometry); CHKERRQ(ierr);
393 }
394 }
395 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
396 for(;mit!=neumann_forces.end();mit++) {
397 ierr = m_field.loop_finite_elements("POTENTIAL_PROBLEM", mit->first, mit->second->getLoopFe()); CHKERRQ(ierr);
398 }
399
400 LaplacianElem elem(m_field,A,F);
401 ierr = m_field.loop_finite_elements("POTENTIAL_PROBLEM", "POTENTIAL_ELEM", elem); CHKERRQ(ierr);
402
403 //post proces fix boundary conditiond
404 ierr = m_field.problem_basic_method_postProcess("POTENTIAL_PROBLEM", fix_dofs); CHKERRQ(ierr);
405
406 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
407 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
408
409 ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
410 ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
411
412 //set matrix possitives define and symetric for cholesky and icc preceonditionser
413 ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE); CHKERRQ(ierr);
414
415 KSP solver;
416 ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
417 ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
418 ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
419 ierr = KSPSetUp(solver); CHKERRQ(ierr);
420
421 ierr = KSPSolve(solver,F,D); CHKERRQ(ierr);
422 ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
423 ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
424 ierr = m_field.set_global_ghost_vector("POTENTIAL_PROBLEM",ROW,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
425
426 if(pcomm->rank()==0) {
427 ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset_fibres); CHKERRQ(ierr);
428 }
429
430 //remove intities from finite elements (to prepare it for next fibres)
431 ierr = m_field.remove_ents_from_finite_element("POTENTIAL_ELEM", new_tets); CHKERRQ(ierr);
432 ierr = m_field.remove_ents_from_finite_element("FLUX_FE", new_tris_1); CHKERRQ(ierr);
433 ierr = m_field.remove_ents_from_finite_element("FLUX_FE", new_tris_2); CHKERRQ(ierr);
434 //clear problems (to prepare it for next fibre)
435 ierr = m_field.clear_problems(); CHKERRQ(ierr);
436
437 ierr = VecDestroy(&F); CHKERRQ(ierr);
438 ierr = VecDestroy(&D); CHKERRQ(ierr);
439 ierr = MatDestroy(&A); CHKERRQ(ierr);
440 ierr = KSPDestroy(&solver); CHKERRQ(ierr);
441 }
442
443
444 Tag th_phi;
445 double def_val = 0;
446 rval = moab.tag_get_handle("PHI",1,MB_TYPE_DOUBLE,th_phi,MB_TAG_CREAT|MB_TAG_SPARSE,&def_val); CHKERRQ_MOAB(rval);
447 for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(m_field,"POTENTIAL_FIELD",dof)) {
448 EntityHandle ent = dof->get()->getEnt();
449 double val = dof->get()->getFieldData();
450 rval = moab.tag_set_data(th_phi,&ent,1,&val); CHKERRQ_MOAB(rval);
451 }
452
453 ProjectionFieldOn10NodeTet ent_method_phi_on_10nodeTet(m_field,"POTENTIAL_FIELD",true,false,"PHI");
454 ierr = m_field.loop_dofs("POTENTIAL_FIELD",ent_method_phi_on_10nodeTet); CHKERRQ(ierr);
455 ent_method_phi_on_10nodeTet.setNodes = false;
456 ierr = m_field.loop_dofs("POTENTIAL_FIELD",ent_method_phi_on_10nodeTet); CHKERRQ(ierr);
457
458 if(pcomm->rank()==0) {
459 rval = moab.write_file("solution_RVE.h5m"); CHKERRQ_MOAB(rval);
460 }
461
462 if(pcomm->rank()==0) {
463 rval = moab.write_file("out_potential_flow.vtk","VTK","",&out_meshset_fibres,1); CHKERRQ_MOAB(rval);
464 }
465
466 }
468
470}
DEPRECATED typedef DirichletFixFieldAtEntitiesBc FixBcAtEntities
DEPRECATED typedef NeumannForcesSurface NeummanForcesSurface
@ COL
@ ROW
#define CATCH_ERRORS
Catch errors.
@ MF_ZERO
Definition definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ H1
continuous field
Definition definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
@ PRESSURESET
@ NODESET
@ SIDESET
@ UNKNOWNNAME
@ BLOCKSET
#define CHKERRQ_MOAB(a)
check error code of MoAB function
constexpr int order
@ F
#define _IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
virtual bool check_field(const std::string &name) const =0
check if field is in database
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.
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
#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.
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode clear_problems(int verb=DEFAULT_VERBOSITY)=0
clear problems
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
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
double D
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
constexpr AssemblyType A
static char help[]
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode remove_ents_from_finite_element(const std::string name, const EntityHandle meshset, const EntityType type, int verb=DEFAULT_VERBOSITY)=0
remove entities from given refinement level to finite element database
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.
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
Deprecated interface functions.
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...
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
Projection of edge entities with one mid-node on hierarchical basis.
virtual MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const =0

Variable Documentation

◆ help

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

Definition at line 44 of file rve_fibre_directions_interface.cpp.