v0.14.0
Loading...
Searching...
No Matches
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
18extern "C" {
19#include <gm_rule.h>
20}
21
22#include <MoFEM.hpp>
23using 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
37using namespace ObosleteUsersModules;
38
39#include <MethodForForceScaling.hpp>
40#include <PotentialFlowFEMethod.hpp>
41#include <SurfacePressure.hpp>
42
43
44
45
46static char help[] = "...\n\n";
47
48int 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}
FixBcAtEntities
DEPRECATED typedef DirichletFixFieldAtEntitiesBc FixBcAtEntities
Definition DirichletBC.hpp:331
NeummanForcesSurface
DEPRECATED typedef NeumannForcesSurface NeummanForcesSurface
Definition SurfacePressure.hpp:1321
main
int main()
Definition adol-c_atom.cpp:46
COL
@ COL
Definition definitions.h:136
ROW
@ ROW
Definition definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MF_ZERO
@ MF_ZERO
Definition definitions.h:111
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
H1
@ H1
continuous field
Definition definitions.h:85
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
PRESSURESET
@ PRESSURESET
Definition definitions.h:165
NODESET
@ NODESET
Definition definitions.h:159
SIDESET
@ SIDESET
Definition definitions.h:160
UNKNOWNNAME
@ UNKNOWNNAME
Definition definitions.h:171
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition definitions.h:467
order
constexpr int order
Definition dg_projection.cpp:18
F
@ F
Definition free_surface.cpp:394
gm_rule.h
_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_
#define _IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(MFIELD, NAME, IT)
Definition Interface.hpp:1878
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::CoreInterface::add_ents_to_field_by_type
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::CoreInterface::problem_basic_method_postProcess
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition MeshsetsManager.hpp:49
MoFEM::CoreInterface::get_problem_finite_elements_entities
virtual MoFEMErrorCode get_problem_finite_elements_entities(const std::string name, const std::string &fe_name, const EntityHandle meshset)=0
add finite elements to the meshset
MoFEM::CoreInterface::add_problem
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
MoFEM::CoreInterface::modify_problem_ref_level_add_bit
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEM::CoreInterface::modify_problem_add_finite_element
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
D
double D
Definition initial_diffusion.cpp:44
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition Exceptions.hpp:74
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
A
constexpr AssemblyType A
Definition operators_tests.cpp:30
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
help
static char help[]
Definition rve_fibre_directions.cpp:46
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::problem_basic_method_preProcess
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
Definition DeprecatedCoreInterface.cpp:18
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
Generated by Doxygen 1.12.0 and hosted at University of Glasgow