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 <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[] )

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 48 of file rve_fibre_directions_interface.cpp.

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

Generated by Doxygen 1.12.0 and hosted at University of Glasgow