v0.13.1
rve_fibre_directions_interface.cpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15//#include "FunctionsForFieldData.hpp"
16//#include "cholesky.hpp"
17
18extern "C" {
19#include <gm_rule.h>
20}
21
22#include <MoFEM.hpp>
23using namespace MoFEM;
24
25
26#include <boost/math/special_functions/round.hpp>
28#include <DirichletBC.hpp>
29
31#include <petsctime.h>
32
33// #include <FEMethod_LowLevelStudent.hpp>
34// #include <FEMethod_UpLevelStudent.hpp>
35// #include <PostProcVertexMethod.hpp>
36// #include <PostProcDisplacementAndStrainOnRefindedMesh.hpp>
37// using namespace ObosleteUsersModules;
38
40// #include "PotentialFlowFEMethod.hpp"
41#include "SurfacePressure.hpp"
42
43
44static char help[] = "...\n\n";
45
46int main(int argc, char *argv[]) {
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
Definition: definitions.h:136
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
@ PRESSURESET
Definition: definitions.h:165
@ NODESET
Definition: definitions.h:159
@ SIDESET
Definition: definitions.h:160
@ UNKNOWNNAME
Definition: definitions.h:171
@ BLOCKSET
Definition: definitions.h:161
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
#define _IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
Definition: Interface.hpp:1833
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 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 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_col(const std::string &fe_name, const std::string &name_row)=0
set field col 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 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
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
char mesh_file_name[255]
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
double A
int main(int argc, char *argv[])
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:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
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