Getting No. of Fibres and their index to be used for Potential Flow Problem
46 {
47
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;
60 if(flg != PETSC_TRUE) {
61 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
62 }
65 if(flg != PETSC_TRUE) {
67 }
68
69
70
71 PetscInt mesh_refinement_level;
73 if(flg != PETSC_TRUE) {
74 mesh_refinement_level = 0;
75 }
76
77
78 const char *option;
79 option = "";
81 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
82 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
83
84
89
90
91
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;
102 cout<<" All nodes on negative surfaces " << SurNodesNeg.size()<<endl;
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++) {
110
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;
115
116 }
117
118 for(Range::iterator nit = SurNodesPos.begin(); nit!=SurNodesPos.end(); nit++) {
120
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;
126
127 }
128
129
130
131
132
133
136
137
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());
158
159 aaa += 1;
160 }
161
162
168
169
171
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];
181
183 vector<BitRefLevel> bit_levels;
185
186
187 int ll = 1;
188
189
190
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);
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
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;
208
209 ierr = interface_ptr->
getSides(cubit_meshset,bit_levels.back(),
true); CHKERRQ(
ierr);
210
212
213 ierr = interface_ptr->
splitSides(ref_level_meshset,bit_levels.back(),cubit_meshset,
true,
true,0); CHKERRQ(
ierr);
214
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
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
231 int meshset_data_root[2]={ll,0};
233
234
237 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),
BitRefLevel().set(),MBTET,out_tet_meshset); CHKERRQ(
ierr);
239
240
241
244 ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),
BitRefLevel().set(),MBPRISM,out_meshset_prism); CHKERRQ(
ierr);
246
247
250 ierr = m_field.get_entities_by_ref_level(bit_levels.back(),
BitRefLevel().set(),out_meshset); CHKERRQ(
ierr);
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
261
262
265
266
271
276
277
284
285
286
290
292 if(pcomm->rank()==0) {
294 }
295
296
297 for (int cc = 0; cc < noOfFibres; cc++) {
298 ostringstream rrr, ppp1, ppp2;
299 rrr << "PotentialFlow_" << fibreList[cc];
302 if(it->getName() == rrr.str().c_str() ) {
303
306 new_tets = intersect(LatestRefinedTets,TetsInBlock);
307
309 }
310 }
311
312 Range new_tris_1, new_tris_2;
313
314 ppp1 << "PressureIO_" << fibreList[cc] <<"_1";
315 ppp2 << "PressureIO_" << fibreList[cc] <<"_2";
317 if (ppp1.str().c_str()==it->getName()){
320 new_tris_1 = intersect(LatestRefinedTris,tris);
322 }
323
324
325 if (ppp2.str().c_str()==it->getName()){
328 new_tris_2 = intersect(LatestRefinedTris,tris);
330 }
331 }
332
333
335
336
338
340
342
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
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
361 if (zero_pressure.str().c_str()==it->getName()){
369 fix_nodes.insert(adj.begin(),adj.end());
371 fix_nodes.insert(adj.begin(),adj.end());
373 }
374 }
376
378
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
385
386 boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
387 string fe_name_str = "FLUX_FE";
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++) {
398 }
399
400 LaplacianElem elem(m_field,A,
F);
402
403
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);
411
412
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
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) {
428 }
429
430
434
436
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);
449 double val = dof->get()->getFieldData();
451 }
452
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) {
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
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERRQ_MOAB(a)
check error code of MoAB function
#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
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.
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)
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
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