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