v0.11.1
Functions | Variables
forces_and_sources_getting_mult_H1_H1_atom_test.cpp File Reference

Atom test verifying forces and sources operator on H1 approx. space. More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions

template<typename T >
void zero_entries (T &t)
 
int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 
constexpr double eps = 1e-6
 

Detailed Description

Atom test verifying forces and sources operator on H1 approx. space.

Definition in file forces_and_sources_getting_mult_H1_H1_atom_test.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 37 of file forces_and_sources_getting_mult_H1_H1_atom_test.cpp.

37  {
38 
39  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
40 
41  try {
42 
43  moab::Core mb_instance;
44  moab::Interface &moab = mb_instance;
45  int rank;
46  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
47 
48  PetscBool flg = PETSC_TRUE;
49  char mesh_file_name[255];
50 #if PETSC_VERSION_GE(3, 6, 4)
51  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
52  255, &flg);
53 #else
54  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
55  mesh_file_name, 255, &flg);
56 #endif
57  if (flg != PETSC_TRUE) {
58  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
59  }
60 
61  // Create MoFEM database
62  MoFEM::Core core(moab);
63  MoFEM::Interface &m_field = core;
64 
65  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
66  if (pcomm == NULL)
67  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
68 
69  const char *option;
70  option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
72  CHKERR moab.load_file(mesh_file_name, 0, option);
73  CHKERRG(rval);
75 
76  // set entitities bit level
77  BitRefLevel bit_level0;
78  bit_level0.set(0);
79  EntityHandle meshset_level0;
80  CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
81  CHKERRG(rval);
82  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
83  0, 3, bit_level0);
84 
85  // Fields
86  CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 1);
87  CHKERR m_field.add_field("FIELD2", H1, AINSWORTH_LEGENDRE_BASE, 3);
88 
89  // FE
90  CHKERR m_field.add_finite_element("TEST_FE");
91 
92  // Define rows/cols and element data
93  CHKERR m_field.modify_finite_element_add_field_row("TEST_FE", "FIELD1");
94  CHKERR m_field.modify_finite_element_add_field_col("TEST_FE", "FIELD2");
95  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE", "FIELD1");
96  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE", "FIELD2");
97 
98  // Problem
99  CHKERR m_field.add_problem("TEST_PROBLEM");
100 
101  // set finite elements for problem
102  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TEST_FE");
103  // set refinement level for problem
104  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
105 
106  // meshset consisting all entities in mesh
107  EntityHandle root_set = moab.get_root_set();
108  // add entities to field
109  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FIELD1");
110  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FIELD2");
111  // add entities to finite element
112  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBTET,
113  "TEST_FE");
114 
115  // set app. order
116  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
117  // (Mark Ainsworth & Joe Coyle)
118  int order = 5;
119  CHKERR m_field.set_field_order(root_set, MBTET, "FIELD1", order);
120  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD1", order);
121  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD1", order);
122  CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD1", 1);
123  CHKERR m_field.set_field_order(root_set, MBTET, "FIELD2", order);
124  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD2", order);
125  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD2", order);
126  CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD2", 1);
127 
128  /****/
129  // build database
130  // build field
131  CHKERR m_field.build_fields();
132  // build finite elemnts
133  CHKERR m_field.build_finite_elements();
134  // build adjacencies
135  CHKERR m_field.build_adjacencies(bit_level0);
136  // build problem
137  //
138  ProblemsManager *prb_mng_ptr;
139  CHKERR m_field.getInterface(prb_mng_ptr);
140  // const Problem_multiIndex *problems_ptr;
141  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
142 
143  /****/
144  // mesh partitioning
145  // partition
146  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
147  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
148  // what are ghost nodes, see Petsc Manual
149  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
150 
151  struct ForcesAndSourcesCore_TestFE : public ForcesAndSourcesCore {
152 
153  typedef tee_device<std::ostream, std::ofstream> TeeDevice;
154  typedef stream<TeeDevice> TeeStream;
155 
156  struct my_mult_H1_H1 : public DataOperator {
157 
158  std::ofstream ofs;
159  TeeDevice my_tee;
160  TeeStream my_split;
161 
162  my_mult_H1_H1()
163  : ofs("forces_and_sources_getting_mult_H1_H1_atom_test.txt"),
164  my_tee(std::cout, ofs), my_split(my_tee){};
165 
166  ~my_mult_H1_H1() { my_split.close(); }
167 
168  ublas::matrix<FieldData> NN;
169 
170  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
171  EntityType col_type,
175 
176  row_data.getBase() = AINSWORTH_LEGENDRE_BASE;
177  col_data.getBase() = AINSWORTH_LEGENDRE_BASE;
178  int nb_row_dofs = row_data.getN().size2();
179  int nb_col_dofs = col_data.getN().size2();
180 
181  my_split << row_side << " " << col_side << " " << row_type << " "
182  << col_type << std::endl;
183  my_split << "nb_row_dofs " << nb_row_dofs << " nb_col_dofs "
184  << nb_col_dofs << std::endl;
185  NN.resize(nb_row_dofs, nb_col_dofs);
186 
187  zero_entries(row_data.getN().data());
188  zero_entries(row_data.getDiffN().data());
189 
190  my_split << std::setprecision(3);
191  my_split << std::fixed;
192  my_split << row_data.getN() << std::endl;
193  my_split << col_data.getN() << std::endl;
194 
195  for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
196 
197  bzero(&*NN.data().begin(),
198  nb_row_dofs * nb_col_dofs * sizeof(FieldData));
199 
200  cblas_dger(CblasRowMajor, nb_row_dofs, nb_col_dofs, 1,
201  &row_data.getN()(gg, 0), 1, &col_data.getN()(gg, 0), 1,
202  &*NN.data().begin(), nb_col_dofs);
203 
204  my_split << "gg " << gg << " : ";
205  my_split << std::setprecision(3);
206  my_split << std::fixed;
207 
208  MatrixDouble difference =
209  NN - outer_prod(row_data.getN(gg), col_data.getN(gg));
210  zero_entries(difference.data());
211 
212  my_split << difference << std::endl;
213  if (row_type != MBVERTEX) {
214  my_split << row_data.getDiffN(gg) << std::endl;
215  }
216 
217  if (row_type == MBVERTEX) {
218  my_split << row_data.getDiffN() << std::endl;
219  } else {
220  typedef ublas::array_adaptor<FieldData> storage_t;
221  storage_t st(nb_row_dofs * 3, &row_data.getDiffN()(gg, 0));
222  ublas::matrix<FieldData, ublas::row_major, storage_t>
223  digNatGaussPt(nb_row_dofs, 3, st);
224  my_split << std::endl << digNatGaussPt << std::endl;
225  }
226  }
227 
228  my_split << std::endl;
229 
231  }
232  };
233 
234  my_mult_H1_H1 op;
235 
236  ForcesAndSourcesCore_TestFE(MoFEM::Interface &_m_field)
237  : ForcesAndSourcesCore(_m_field), data_row(MBTET), data_col(MBTET){};
238 
242  }
243 
244  DataForcesAndSourcesCore data_row, data_col;
245 
248 
251 
252  CHKERR getEntitySense<MBEDGE>(data_row);
253  CHKERR getEntitySense<MBTRI>(data_row);
254  CHKERR getEntitySense<MBEDGE>(data_col);
255  CHKERR getEntitySense<MBTRI>(data_col);
256 
257  CHKERR getEntityDataOrder<MBEDGE>(data_row, H1);
258  CHKERR getEntityDataOrder<MBEDGE>(data_col, H1);
259  CHKERR getEntityDataOrder<MBTRI>(data_row, H1);
260  CHKERR getEntityDataOrder<MBTRI>(data_col, H1);
261  CHKERR getEntityDataOrder<MBTET>(data_row, H1);
262  CHKERR getEntityDataOrder<MBTET>(data_col, H1);
263  data_row.dataOnEntities[MBVERTEX][0].getBase() =
265  CHKERR getEntityFieldData(data_row, "FIELD1", MBEDGE);
266  data_col.dataOnEntities[MBVERTEX][0].getBase() =
268  CHKERR getEntityFieldData(data_col, "FIELD2", MBEDGE);
269  CHKERR getRowNodesIndices(data_row, "FIELD1");
270  CHKERR getColNodesIndices(data_col, "FIELD2");
271  CHKERR getEntityRowIndices(data_row, "FIELD1", MBEDGE);
272  CHKERR getEntityColIndices(data_col, "FIELD2", MBEDGE);
273  CHKERR getFaceTriNodes(data_row);
274  CHKERR getFaceTriNodes(data_col);
275 
276  MatrixDouble gauss_pts(4, 4);
277  for (int gg = 0; gg < 4; gg++) {
278  gauss_pts(0, gg) = G_TET_X4[gg];
279  gauss_pts(1, gg) = G_TET_Y4[gg];
280  gauss_pts(2, gg) = G_TET_Z4[gg];
281  gauss_pts(3, gg) = G_TET_W4[gg];
282  }
284  gauss_pts,
285  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
286  data_row, H1, AINSWORTH_LEGENDRE_BASE)));
288  gauss_pts,
289  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
290  data_col, H1, AINSWORTH_LEGENDRE_BASE)));
291 
292  CHKERR op.opLhs(data_row, data_col);
293 
295  }
296 
299 
301  }
302  };
303 
304  ForcesAndSourcesCore_TestFE fe1(m_field);
305  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TEST_FE", fe1);
306  }
307  CATCH_ERRORS;
308 
310 
311  return 0;
312 }
@ CblasRowMajor
Definition: cblas.h:10
void cblas_dger(const enum CBLAS_ORDER order, const int M, const int N, const double alpha, const double *X, const int incX, const double *Y, const int incY, double *A, const int lda)
Definition: cblas_dger.c:12
#define BARRIER_PCOMM_RANK_END(PCMB)
set barrier start Run code in sequence, starting from process 0, and ends on last process.
Definition: definitions.h:342
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
@ H1
continuous field
Definition: definitions.h:177
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
#define BARRIER_PCOMM_RANK_START(PCMB)
set barrier start Run code in sequence, starting from process 0, and ends on last process.
Definition: definitions.h:320
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
static const double G_TET_W4[]
Definition: fem_tools.h:1098
static const double G_TET_Z4[]
Definition: fem_tools.h:1096
static const double G_TET_Y4[]
Definition: fem_tools.h:1094
static const double G_TET_X4[]
Definition: fem_tools.h:1092
tee_device< std::ostream, std::ofstream > TeeDevice
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
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_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
double FieldData
Field data type.
Definition: Types.hpp:36
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
CoreTmp< 0 > Core
Definition: Core.hpp:1140
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1953
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
FieldApproximationBase & getBase()
Get approximation base.
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
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)
Class used to pass element data to calculate base functions on tet,triangle,edge.
structure to get information form mofem into DataForcesAndSourcesCore
MoFEMErrorCode getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getColNodesIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get col node indices from FENumeredDofEntity_multiIndex
MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const
Get nodes on triangles.
virtual MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode getEntityRowIndices(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode getEntityColIndices(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode getRowNodesIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get row node indices from FENumeredDofEntity_multiIndex
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
Problem manager is used to build and partition problems.
Calculate base functions on tetrahedral.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.

◆ zero_entries()

template<typename T >
void zero_entries ( T t)

Definition at line 31 of file forces_and_sources_getting_mult_H1_H1_atom_test.cpp.

31  {
32  for (auto &v : t)
33  if (std::abs(v) < eps)
34  v = 0;
35 }

Variable Documentation

◆ eps

constexpr double eps = 1e-6
constexpr

◆ help

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