v0.10.0
Functions | Variables
best_approximation_obsolete.cpp File Reference
#include <MoFEM.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <FieldApproximation.hpp>
#include <PotsProcOnRefMesh.hpp>
#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>
#include <petsctime.h>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <cmath>
#include <boost/math/special_functions.hpp>
#include <complex>
#include <AnalyticalSolutions.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[] 
)

Definition at line 52 of file best_approximation_obsolete.cpp.

52  {
53 
54  ErrorCode rval;
55  PetscErrorCode ierr;
56 
57  PetscInitialize(&argc,&argv,(char *)0,help);
58 
59  moab::Core mb_instance;
60  Interface& moab = mb_instance;
61  int rank;
62  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
63 
64  PetscBool flg = PETSC_TRUE;
65  char mesh_file_name[255];
66  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
67  if(flg != PETSC_TRUE) {
68  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
69  }
70 
71  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
72  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
73 
74  PetscBool is_partitioned = PETSC_FALSE;
75  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
76  if(is_partitioned == PETSC_TRUE) {
77  //Read mesh to MOAB
78  const char *option;
79  option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
80  rval = moab.load_file(mesh_file_name, 0, option); CHKERR_PETSC(rval);
81  } else {
82  const char *option;
83  option = "";
84  rval = moab.load_file(mesh_file_name, 0, option); CHKERR_PETSC(rval);
85  }
86 
87  // create MoFEM database
88  MoFEM::Core core(moab);
89  MoFEM::Interface& m_field = core;
90 
91  // count the comsumption of time by single run
92  PetscLogDouble t1,t2;
93  PetscLogDouble v1,v2;
94  ierr = PetscTime(&v1); CHKERRQ(ierr);
95  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
96 
97  // set entitities bit level
98  BitRefLevel bit_level0;
99  bit_level0.set(0);
100  EntityHandle meshset_level0;
101  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERR_PETSC(rval);
102  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
103 
104  // define fields
105  ierr = m_field.add_field("reEX",H1,1); CHKERRQ(ierr);
106  ierr = m_field.add_field("imEX",H1,1); CHKERRQ(ierr);
107  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,3,MF_ZERO); CHKERRQ(ierr);
108 
109  // define finite element
110  ierr = m_field.add_finite_element("FE1"); CHKERRQ(ierr);
111 
112  // Define rows/cols and element data
113  ierr = m_field.modify_finite_element_add_field_row("FE1","reEX"); CHKERRQ(ierr);
114  ierr = m_field.modify_finite_element_add_field_col("FE1","reEX"); CHKERRQ(ierr);
115  ierr = m_field.modify_finite_element_add_field_data("FE1","reEX"); CHKERRQ(ierr);
116  ierr = m_field.modify_finite_element_add_field_data("FE1","imEX"); CHKERRQ(ierr);
117  ierr = m_field.modify_finite_element_add_field_data("FE1","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
118 
119  if(m_field.check_field("rePRES") && m_field.check_field("imPRES")) {
120 
121  ierr = m_field.modify_finite_element_add_field_data("FE1","rePRES"); CHKERRQ(ierr);
122  ierr = m_field.modify_finite_element_add_field_data("FE1","imPRES"); CHKERRQ(ierr);
123 
124  }
125 
126  // meshset consisting all entities in mesh
127  EntityHandle root_set = moab.get_root_set();
128  // add entities to field
129  ierr = m_field.add_ents_to_field_by_TETs(root_set,"reEX"); CHKERRQ(ierr);
130  ierr = m_field.add_ents_to_field_by_TETs(root_set,"imEX"); CHKERRQ(ierr);
131  ierr = m_field.add_ents_to_field_by_TETs(root_set,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
132  // add entities to finite element
133  ierr = m_field.add_ents_to_finite_element_by_TETs(root_set,"FE1"); CHKERRQ(ierr);
134 
135  // set app. order
136  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
137  int order = 3;
138  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
139  if(flg != PETSC_TRUE) {
140  order = 3;
141  }
142  ierr = m_field.set_field_order(root_set,MBTET,"reEX",order); CHKERRQ(ierr);
143  ierr = m_field.set_field_order(root_set,MBTRI,"reEX",order); CHKERRQ(ierr);
144  ierr = m_field.set_field_order(root_set,MBEDGE,"reEX",order); CHKERRQ(ierr);
145  ierr = m_field.set_field_order(root_set,MBVERTEX,"reEX",1); CHKERRQ(ierr);
146 
147  ierr = m_field.set_field_order(root_set,MBTET,"imEX",order); CHKERRQ(ierr);
148  ierr = m_field.set_field_order(root_set,MBTRI,"imEX",order); CHKERRQ(ierr);
149  ierr = m_field.set_field_order(root_set,MBEDGE,"imEX",order); CHKERRQ(ierr);
150  ierr = m_field.set_field_order(root_set,MBVERTEX,"imEX",1); CHKERRQ(ierr);
151 
152  ierr = m_field.set_field_order(root_set,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
153  ierr = m_field.set_field_order(root_set,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
154  ierr = m_field.set_field_order(root_set,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
155  ierr = m_field.set_field_order(root_set,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
156 
157  // build field
158  ierr = m_field.build_fields(); CHKERRQ(ierr);
159  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
160  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
161  // build finite elemnts
162  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
163  // build adjacencies
164  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
165 
166  // define problem
167  ierr = m_field.add_problem("EX1_PROBLEM"); CHKERRQ(ierr);
168  // set finite elements for problem
169  ierr = m_field.modify_problem_add_finite_element("EX1_PROBLEM","FE1"); CHKERRQ(ierr);
170  // set refinment level for problem
171  ierr = m_field.modify_problem_ref_level_add_bit("EX1_PROBLEM",bit_level0); CHKERRQ(ierr);
172 
173  // build porblems
174  if(is_partitioned) {
175  // if mesh is partitioned
176  ierr = m_field.build_partitioned_problem("EX1_PROBLEM",true); CHKERRQ(ierr);
177  ierr = m_field.partition_finite_elements("EX1_PROBLEM",true); CHKERRQ(ierr);
178  } else {
179  // if not partitioned mesh is load to all processes
180  ierr = m_field.build_problem("EX1_PROBLEM"); CHKERRQ(ierr);
181  ierr = m_field.partition_problem("EX1_PROBLEM"); CHKERRQ(ierr);
182  ierr = m_field.partition_finite_elements("EX1_PROBLEM"); CHKERRQ(ierr);
183  }
184  ierr = m_field.partition_ghost_dofs("EX1_PROBLEM"); CHKERRQ(ierr);
185 
186 
187 
188  // extract data from MAT_HELMHOLTZ block
189  double angularfreq = 1;
190  double speed = 1;
191 
192  /// this works only for one block
193  int nb_of_blocks = 0;
194  for(_IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(m_field,"MAT_HELMHOLTZ",it)) {
195 
196  // get block attributes
197  vector<double> attributes;
198  ierr = it->get_Cubit_attributes(attributes); CHKERRQ(ierr);
199  if(attributes.size()<2) {
200  SETERRQ1(PETSC_COMM_SELF,MOFEM_INVALID_DATA,
201  "not enough block attributes, expected 2 attributes ( angular freq., speed) , attributes.size() = %d ",attributes.size());
202  }
203  angularfreq = attributes[0];
204  speed = attributes[1];
205  nb_of_blocks++;
206 
207  }
208 
209  if(nb_of_blocks!=1) {
210  PetscPrintf(PETSC_COMM_SELF,"Warrning: wave number is set to all blocks baesd on last evaluated block");
211  }
212  double wavenumber = angularfreq/speed;
213 
214  // set wave number from line command, that overwrite numbre form block set
215  ierr = PetscOptionsGetScalar(NULL,"-wave_number",&wavenumber,NULL); CHKERRQ(ierr);
216 
217  //wave direction unit vector=[x,y,z]^T
218  double waveDirection[3];
219  int nmax=3;
220  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-wave_direction",waveDirection,&nmax,&flg); CHKERRQ(ierr);
221 
222  VectorDouble wave_direction;
223  wave_direction.resize(3);
224  cblas_dcopy(3, &waveDirection[0], 1, &wave_direction(0), 1);
225  if(flg != PETSC_TRUE) {
226  wave_direction[0] = 1;
227  wave_direction[1] = 0;
228  wave_direction[2] = 0;
229  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -wave_direction [3*1 vector] default:X direction [1,0,0]");
230  }
231 
232  PetscInt choise_value = 0;
233  // set type of analytical solution
234  ierr = PetscOptionsGetEList(NULL,"-analytical_solution_type",analytical_solution_types,6,&choise_value,&flg); CHKERRQ(ierr);
235  if(flg != PETSC_TRUE) {
236  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -analytical_solution_type needed, WARNING!!!!!!.");
237  }
238 
239 
240  switch((AnalyticalSolutionTypes)choise_value) {
241 
243 
244  {
245  HardSphereScatterWave function_evaluator(wavenumber);
246  ierr = solve_problem(m_field,"EX1_PROBLEM","FE1","reEX","imEX",INSERT_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
247  }
248 
249  break;
250 
251 
253 
254  {
255  SoftSphereScatterWave function_evaluator(wavenumber);
256  ierr = solve_problem(m_field,"EX1_PROBLEM","FE1","reEX","imEX",INSERT_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
257  }
258 
259  break;
260 
261  case PLANE_WAVE:
262 
263  {
264  PlaneWave function_evaluator(wavenumber,0.25*M_PI);
265  ierr = solve_problem(m_field,"EX1_PROBLEM","FE1","reEX","imEX",INSERT_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
266  }
267 
268  break;
269 
271 
272  {
273  HardCylinderScatterWave function_evaluator(wavenumber);
274  ierr = solve_problem(m_field,"EX1_PROBLEM","FE1","reEX","imEX",INSERT_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
275  }
276 
277  break;
278 
280 
281  {
282  SoftCylinderScatterWave function_evaluator(wavenumber);
283  ierr = solve_problem(m_field,"EX1_PROBLEM","FE1","reEX","imEX",INSERT_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
284  }
285 
286  break;
287 
288  case INCIDENT_WAVE:
289 
290  {
291  IncidentWave function_evaluator(wavenumber,wave_direction);
292  ierr = solve_problem(m_field,"EX1_PROBLEM","FE1","reEX","imEX",INSERT_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
293  }
294 
295  break;
296 
297  }
298 
299  PetscBool add_incident_wave = PETSC_FALSE;
300  ierr = PetscOptionsGetBool(NULL,"-add_incident_wave",&add_incident_wave,NULL); CHKERRQ(ierr);
301  if(add_incident_wave) {
302 
303  IncidentWave function_evaluator(wavenumber,wave_direction);
304  ierr = solve_problem(m_field,"EX1_PROBLEM","FE1","reEX","imEX",ADD_VALUES,function_evaluator,is_partitioned); CHKERRQ(ierr);
305 
306  }
307 
308  if(is_partitioned) {
309  rval = moab.write_file("best_solution.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
310  } else {
311  if(!pcomm->rank()) {
312  rval = moab.write_file("analytical_solution.h5m"); CHKERR_PETSC(rval);
313  }
314  }
315 
316 
317  PetscBool save_postproc_mesh = PETSC_FALSE;
318  ierr = PetscOptionsGetBool(NULL,"-save_postproc_mesh",&save_postproc_mesh,NULL); CHKERRQ(ierr);
319  if(save_postproc_mesh) {
320 
322  ierr = post_proc.generateRefereneElemenMesh(); CHKERRQ(ierr);
323  ierr = post_proc.addFieldValuesPostProc("reEX"); CHKERRQ(ierr);
324  ierr = post_proc.addFieldValuesPostProc("imEX"); CHKERRQ(ierr);
325 
326  if(m_field.check_field("rePRES") && m_field.check_field("imPRES")) {
327  ierr = post_proc.addFieldValuesPostProc("rePRES"); CHKERRQ(ierr);
328  ierr = post_proc.addFieldValuesPostProc("imPRES"); CHKERRQ(ierr);
329  }
330 
331  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
332  ierr = m_field.loop_finite_elements("EX1_PROBLEM","FE1",post_proc); CHKERRQ(ierr);
333  rval = post_proc.postProcMesh.write_file("best_solution_mesh_post_proc.h5m","MOAB","PARALLEL=WRITE_PART"); CHKERR_PETSC(rval);
334 
335  }
336 
337  // calulate total time
338  ierr = PetscTime(&v2);CHKERRQ(ierr);
339  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
340  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f S CPU Time = %f S \n",pcomm->rank(),v2-v1,t2-t1);
341 
342  ierr = PetscFinalize(); CHKERRQ(ierr);
343 
344  return 0;
345 
346 }
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
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
CoreTmp< 0 > Core
Definition: Core.hpp:1129
Calculate the analytical solution of impinging wave on cylinder.
auto post_proc
static char help[]
#define _IT_CUBITMESHSETS_BY_NAME_FOR_LOOP_(MESHSET_MANAGER, NAME, IT)
Iterator that loops over Cubit BlockSet having a particular name.
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
virtual MoFEMErrorCode 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.
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.
Calculate the analytical solution of impinging wave on cylinder.
char mesh_file_name[255]
Projection of edge entities with one mid-node on hierarchical basis.
Calculate the analytical solution of plane wave guide propagating in direction theta.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
PetscErrorCode solve_problem(MoFEM::Interface &m_field, const string &problem_name, const string &fe_name, const string &re_field, const string &im_field, InsertMode mode, FUNEVAL &fun_evaluator, PetscBool is_partitioned)
const char * analytical_solution_types[]
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DEPRECATED MoFEMErrorCode partition_finite_elements(const std::string &name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=-1)
partition finite elements
constexpr int order
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
CHKERRQ(ierr)
DEPRECATED MoFEMErrorCode build_problem(const std::string &name, const bool square_matrix, int verb=-1)
build problem data structures
Post processing.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
DEPRECATED MoFEMErrorCode partition_problem(const std::string &name, int verb=-1)
partition problem dofs (collective)
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
AnalyticalSolutionTypes
virtual DEPRECATED MoFEMErrorCode add_ents_to_finite_element_by_TETs(const Range &tets, const std::string &name)=0
add TET entities from range to finite element database given by name
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
DEPRECATED MoFEMErrorCode partition_ghost_dofs(const std::string &name, int verb=-1)
determine ghost nodes
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual bool check_field(const std::string &name) const =0
check if field is in database
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
Core (interface) class.
Definition: Core.hpp:77
continuous field
Definition: definitions.h:177
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
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...
DEPRECATED MoFEMErrorCode add_ents_to_field_by_TETs(const EntityHandle meshset, const std::string &name, int verb=-1)
set field entities from adjacencies of tetrahedronThe lower dimension entities are added depending on...
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)

Variable Documentation

◆ help

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

Definition at line 49 of file best_approximation_obsolete.cpp.