v0.14.0
Functions | Variables
rve_fibre_directions.cpp File Reference
#include <gm_rule.h>
#include <MoFEM.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[] 
)

Populating the Multi-index container with nodes on +ve faces

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

Adding entities to Field and FE for Potential Flow Problem

Definition at line 48 of file rve_fibre_directions.cpp.

48  {
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 }

Variable Documentation

◆ help

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

Definition at line 46 of file rve_fibre_directions.cpp.

SIDESET
@ SIDESET
Definition: definitions.h:147
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:215
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:454
PRESSURESET
@ PRESSURESET
Definition: definitions.h:152
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
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:158
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
NODESET
@ NODESET
Definition: definitions.h:146
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
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
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.
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.
double
COL
@ COL
Definition: definitions.h:123
NeummanForcesSurface
DEPRECATED typedef NeumannForcesSurface NeummanForcesSurface
Definition: SurfacePressure.hpp:1321
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
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:23
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:98
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_
#define _IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
Definition: Interface.hpp:1843
help
static char help[]
Definition: rve_fibre_directions.cpp:46
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
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
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_filed)=0
set finite element field data
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.
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
F
@ F
Definition: free_surface.cpp:394
Generated by Doxygen 1.8.17 and hosted at University of Glasgow