v0.9.0
Functions | Variables
check_base_functions_direvatives_tri.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Functions

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

Variables

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

Function Documentation

◆ main()

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

Definition at line 28 of file check_base_functions_direvatives_tri.cpp.

28  {
29 
30  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
31 
32  try {
33 
34  enum spaces { H1TRI, HCURLTRI, LASTOP };
35 
36  const char *list_spaces[] = {"h1tri", "hcurltri"};
37 
38  PetscBool flg;
39  PetscInt choice_value = H1TRI;
40  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces, LASTOP,
41  &choice_value, &flg);
42 
43  if (flg != PETSC_TRUE) {
44  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
45  }
46 
47  FieldSpace space = LASTSPACE;
48  if (choice_value == H1TRI) {
49  space = H1;
50  } else if (choice_value == HCURLTRI) {
51  space = HCURL;
52  }
53 
54  // Select base
55  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
56 
57  const char *list_bases[] = {"ainsworth", "demkowicz"};
58 
59  PetscInt choice_base_value = AINSWORTH;
60  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
61  LASBASETOP, &choice_base_value, &flg);
62 
63  if (flg != PETSC_TRUE) {
64  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
65  }
66 
68  if (choice_base_value == AINSWORTH) {
70  } else if (choice_base_value == DEMKOWICZ) {
71  base = DEMKOWICZ_JACOBI_BASE;
72  }
73 
74 
75  moab::Core mb_instance;
76  moab::Interface &moab = mb_instance;
77  int rank;
78  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
79 
80  // create one tet
81  double tri_coords[] = {0, 0, 0,
82 
83  0.5, 0, 0,
84 
85  0, 2., 0};
86  EntityHandle nodes[3];
87  for (int nn = 0; nn < 3; nn++) {
88  CHKERR moab.create_vertex(&tri_coords[3 * nn], nodes[nn]);
89 
90  }
91  EntityHandle tri;
92  CHKERR moab.create_element(MBTRI, nodes, 3, tri);
93 
94  // Create adjacencies entities
95  Range adj;
96  CHKERR moab.get_adjacencies(&tri, 1, 1, true, adj);
97 
98 
99  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
100  if (pcomm == NULL)
101  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
102 
103  // Create MoFEM database
104  MoFEM::Core core(moab);
105  MoFEM::Interface &m_field = core;
106 
107  // set entities bit level
108  BitRefLevel bit_level0;
109  bit_level0.set(0);
110  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
111  0, 2, bit_level0);
112 
113 
114  // Fields
115  CHKERR m_field.add_field("FIELD", space, base, 1);
116 
117 
118  // FE TET
119  CHKERR m_field.add_finite_element("TRI_FE");
120 
121  // Define rows/cols and element data
122  CHKERR m_field.modify_finite_element_add_field_row("TRI_FE", "FIELD");
123  CHKERR m_field.modify_finite_element_add_field_col("TRI_FE", "FIELD");
124  CHKERR m_field.modify_finite_element_add_field_data("TRI_FE", "FIELD");
125 
126 
127  // Problem
128  CHKERR m_field.add_problem("TEST_PROBLEM");
129 
130  // set finite elements for problem
131  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TRI_FE");
132 
133  // set refinement level for problem
134  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
135 
136 
137  // meshset consisting all entities in mesh
138  EntityHandle root_set = moab.get_root_set();
139  // add entities to field
140  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTRI, "FIELD");
141 
142  // add entities to finite element
143  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBTRI,
144  "TRI_FE");
145 
146  // set app. order
147  int order = 5;
148  if (space == H1) {
149  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD", order);
150  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD", order);
151  CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD", 1);
152 
153  }
154  if (space == HCURL) {
155  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD", order);
156  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD", order);
157  }
158 
159  // build field
160  CHKERR m_field.build_fields();
161 
162  // build finite elemnts
163  CHKERR m_field.build_finite_elements();
164 
165  // build adjacencies
166  CHKERR m_field.build_adjacencies(bit_level0);
167 
168  // build problem
169  ProblemsManager *prb_mng_ptr;
170  CHKERR m_field.getInterface(prb_mng_ptr);
171 
172  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
173 
174 
175  // partition
176  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
177  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
178 
179  // what are ghost nodes, see Petsc Manual
180  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
181 
182  typedef tee_device<std::ostream, std::ofstream> TeeDevice;
183  typedef stream<TeeDevice> TeeStream;
184 
185  std::ofstream ofs("forces_and_sources_checking_direvatives.txt");
186  TeeDevice my_tee(std::cout, ofs);
187  TeeStream my_split(my_tee);
188 
189  struct OpCheckingDirevatives
191 
192  TeeStream &mySplit;
193  OpCheckingDirevatives(TeeStream &my_split)
195  "FIELD", UserDataOperator::OPROW),
196  mySplit(my_split) {}
197 
198  MoFEMErrorCode doWork(int side, EntityType type,
201 
202  if (data.getFieldData().size() == 0)
204  mySplit << "type " << type << " side " << side << endl;
205 
206  if (data.getFieldDofs()[0]->getSpace() == H1) {
207 
208  // mySplit << std::fixed << data.getN() << std::endl;
209  // mySplit << std::fixed << data.getDiffN() << std::endl;
210 
211  const int nb_dofs = data.getN().size2();
212  for (int dd = 0; dd != nb_dofs; dd++) {
213  const double dksi =
214  (data.getN()(1, dd) - data.getN()(0, dd)) / (eps);
215  const double deta =
216  (data.getN()(3, dd) - data.getN()(2, dd)) / (4 * eps);
217  mySplit << "DKsi " << dksi << std::endl;
218  mySplit << "DEta " << deta << std::endl;
219  mySplit << "diffN " << data.getDiffN()(4, 2 * dd + 0) << std::endl;
220  mySplit << "diffN " << data.getDiffN()(4, 2 * dd + 1) << std::endl;
221  if (fabs(dksi - data.getDiffN()(4, 2 * dd + 0)) > eps_diff) {
222  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
223  "H1 inconsistent dKsi derivative");
224  }
225  if (fabs(deta - data.getDiffN()(4, 2 * dd + 1)) > eps_diff) {
226  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
227  "H1 inconsistent dEta derivative");
228  }
229  }
230  }
231 
232  if (data.getFieldDofs()[0]->getSpace() == HCURL) {
233 
234  FTensor::Tensor1<double *, 3> base_ksi_m(&data.getN()(0, HVEC0),
235  &data.getN()(0, HVEC1),
236  &data.getN()(0, HVEC2), 3);
237  FTensor::Tensor1<double *, 3> base_ksi_p(&data.getN()(1, HVEC0),
238  &data.getN()(1, HVEC1),
239  &data.getN()(1, HVEC2), 3);
240  FTensor::Tensor1<double *, 3> base_eta_m(&data.getN()(2, HVEC0),
241  &data.getN()(2, HVEC1),
242  &data.getN()(2, HVEC2), 3);
243  FTensor::Tensor1<double *, 3> base_eta_p(&data.getN()(3, HVEC0),
244  &data.getN()(3, HVEC1),
245  &data.getN()(3, HVEC2), 3);
246 
247  // cerr << data.getN() << endl;
248  // cerr << data.getDiffN() << endl;
249 
251  &data.getDiffN()(4, HVEC0_0), &data.getDiffN()(4, HVEC0_1),
252  &data.getDiffN()(4, HVEC1_0), &data.getDiffN()(4, HVEC1_1),
253  &data.getDiffN()(4, HVEC2_0), &data.getDiffN()(4, HVEC2_1));
254 
258 
259  const int nb_dofs = data.getN().size2() / 3;
260  for (int dd = 0; dd != nb_dofs; dd++) {
261 
262  mySplit << "MoFEM " << diff_base(0, 0) << " " << diff_base(1, 0)
263  << " " << diff_base(2, 0) << endl;
264  mySplit << "MoFEM " << diff_base(0, 1) << " " << diff_base(1, 1)
265  << " " << diff_base(2, 1) << endl;
266 
268  dksi(i) = (base_ksi_p(i) - base_ksi_m(i)) / (eps);
270  deta(i) = (base_eta_p(i) - base_eta_m(i)) / (4 * eps);
271  mySplit << "Finite difference dKsi " << dksi(0) << " " << dksi(1)
272  << " " << dksi(2) << endl;
273  mySplit << "Finite difference dEta " << deta(0) << " " << deta(1)
274  << " " << deta(2) << endl;
275 
276  dksi(i) -= diff_base(i, N0);
277  deta(i) -= diff_base(i, N1);
278  if (sqrt(dksi(i) * dksi(i)) > eps_diff) {
279  // mySplit << "KSI ERROR\n";
280  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
281  "%s inconsistent dKsi derivative for type %d",
282  FieldSpaceNames[data.getFieldDofs()[0]->getSpace()],
283  type);
284  } else {
285  mySplit << "OK" << std::endl;
286  }
287  if (sqrt(deta(i) * deta(i)) > eps_diff) {
288  // mySplit << "ETA ERROR\n";
289  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
290  "%s inconsistent dEta derivative for type %d",
291  FieldSpaceNames[data.getFieldDofs()[0]->getSpace()],
292  type);
293  } else {
294  mySplit << "OK" << std::endl;
295  }
296 
297  ++base_ksi_m;
298  ++base_ksi_p;
299  ++base_eta_m;
300  ++base_eta_p;
301  ++diff_base;
302  }
303  }
304 
305  mySplit << endl;
306 
308  }
309  };
310 
311  struct MyFE : public FaceElementForcesAndSourcesCore {
312 
313  MyFE(MoFEM::Interface &m_field)
314  : FaceElementForcesAndSourcesCore(m_field) {}
315  int getRule(int order) { return -1; };
316 
317  MoFEMErrorCode setGaussPts(int order) {
319 
320  const double ksi = G_TRI_X1[0];
321  const double eta = G_TRI_Y1[0];
322 
323  gaussPts.resize(3, 5);
324  gaussPts.clear();
325 
326  gaussPts(0, 0) = ksi - eps;
327  gaussPts(1, 0) = eta;
328  gaussPts(0, 1) = ksi + eps;
329  gaussPts(1, 1) = eta;
330 
331  gaussPts(0, 2) = ksi;
332  gaussPts(1, 2) = eta - eps;
333  gaussPts(0, 3) = ksi;
334  gaussPts(1, 3) = eta + eps;
335 
336  gaussPts(0, 4) = ksi;
337  gaussPts(1, 4) = eta;
338 
339  for (unsigned int ii = 0; ii != gaussPts.size2(); ii++) {
340  gaussPts(2, ii) = 1;
341  }
342 
343  // cerr << gaussPts << endl;
344 
346  }
347  };
348 
349  MyFE tri_fe(m_field);
350 
351  MatrixDouble inv_jac;
352  tri_fe.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
353  if (space == H1) {
354  tri_fe.getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac));
355  }
356  if (space == HCURL) {
357  tri_fe.getOpPtrVector().push_back(new OpSetInvJacHcurlFace(inv_jac));
358  }
359  tri_fe.getOpPtrVector().push_back(new OpCheckingDirevatives(my_split));
360  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TRI_FE", tri_fe);
361 
362  cerr << inv_jac << endl;
363  }
364  CATCH_ERRORS;
365 
367 
368 }
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 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
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
static const double eps_diff
static const double eps
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
static const double G_TRI_X1[]
Definition: fem_tools.h:247
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
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
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
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.
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:175
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Calculate inverse of jacobian for face element.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
FieldApproximationBase
approximation base
Definition: definitions.h:144
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
static const char *const FieldSpaceNames[]
Definition: definitions.h:178
field with continuous tangents
Definition: definitions.h:172
FieldSpace
approximation spaces
Definition: definitions.h:168
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:285
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
brief Transform local reference derivatives of shape function to global derivatives for face
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elementsFunction which partition finite elements based on dofs partitioning....
static const double G_TRI_Y1[]
Definition: fem_tools.h:250
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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
continuous field
Definition: definitions.h:171
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
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...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
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)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions

Variable Documentation

◆ eps

const double eps = 1e-6
static

Definition at line 25 of file check_base_functions_direvatives_tri.cpp.

◆ eps_diff

const double eps_diff = 1e-4
static

Definition at line 26 of file check_base_functions_direvatives_tri.cpp.

◆ help

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

Definition at line 23 of file check_base_functions_direvatives_tri.cpp.