v0.9.0
check_base_functions_direvatives_on_tet.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 <MoFEM.hpp>
16 
17 namespace bio = boost::iostreams;
18 using bio::stream;
19 using bio::tee_device;
20 
21 using namespace MoFEM;
22 
23 static char help[] = "...\n\n";
24 
25 static const double eps = 1e-6;
26 static const double eps_diff = 1e-8;
27 
28 int main(int argc, char *argv[]) {
29 
30  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
31 
32  try {
33 
34  // Select space
35  enum spaces { H1TET, HDIVTET, HCURLTET, LASTSPACEOP };
36 
37  const char *list_spaces[] = {"h1tet", "hdivtet", "hcurltet"};
38 
39  PetscBool flg;
40  PetscInt choise_space_value = H1TET;
41  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
42  LASTSPACEOP, &choise_space_value, &flg);
43 
44  if (flg != PETSC_TRUE) {
45  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "space not set");
46  }
47 
48  FieldSpace space = LASTSPACE;
49  if (choise_space_value == H1TET) {
50  space = H1;
51  } else if (choise_space_value == HDIVTET) {
52  space = HDIV;
53  } else if (choise_space_value == HCURLTET) {
54  space = HCURL;
55  }
56 
57  // Select base
58  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
59 
60  const char *list_bases[] = {"ainsworth", "demkowicz"};
61 
62  PetscInt choice_base_value = AINSWORTH;
63  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
64  LASBASETOP, &choice_base_value, &flg);
65 
66  if (flg != PETSC_TRUE) {
67  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
68  }
69 
71  if (choice_base_value == AINSWORTH) {
73  } else if (choice_base_value == DEMKOWICZ) {
74  base = DEMKOWICZ_JACOBI_BASE;
75  }
76 
77  moab::Core mb_instance;
78  moab::Interface &moab = mb_instance;
79  int rank;
80  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
81 
82  // create one tet
83  double tet_coords[] = {0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5};
84  EntityHandle nodes[4];
85  for (int nn = 0; nn < 4; nn++) {
86  CHKERR moab.create_vertex(&tet_coords[3 * nn], nodes[nn]);
87  CHKERRG(rval);
88  }
89  EntityHandle tet;
90  CHKERR moab.create_element(MBTET, nodes, 4, tet);
91  CHKERRG(rval);
92 
93  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
94  if (pcomm == NULL)
95  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
96 
97  // Create MoFEM database
98  MoFEM::Core core(moab);
99  MoFEM::Interface &m_field = core;
100 
101  // set entities bit level
102  BitRefLevel bit_level0;
103  bit_level0.set(0);
104  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
105  0, 3, bit_level0);
106 
107 
108  // Fields
109  CHKERR m_field.add_field("FIELD", space, base, 1);
110 
111  // FE TET
112  CHKERR m_field.add_finite_element("TET_FE");
113 
114  // Define rows/cols and element data
115  CHKERR m_field.modify_finite_element_add_field_row("TET_FE", "FIELD");
116  CHKERR m_field.modify_finite_element_add_field_col("TET_FE", "FIELD");
117  CHKERR m_field.modify_finite_element_add_field_data("TET_FE", "FIELD");
118 
119 
120  // Problem
121  CHKERR m_field.add_problem("TEST_PROBLEM");
122 
123 
124  // set finite elements for problem
125  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TET_FE");
126 
127  // set refinement level for problem
128  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
129 
130 
131  // meshset consisting all entities in mesh
132  EntityHandle root_set = moab.get_root_set();
133  // add entities to field
134  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FIELD");
135 
136  // add entities to finite element
137  CHKERR
138  m_field.add_ents_to_finite_element_by_type(root_set, MBTET, "TET_FE");
139 
140 
141  // set app. order
142  int order = 5;
143  if (space == H1) {
144  CHKERR m_field.set_field_order(root_set, MBTET, "FIELD", order);
145  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD", order);
146  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD", order);
147  CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD", 1);
148 
149  }
150  if (space == HCURL) {
151  CHKERR m_field.set_field_order(root_set, MBTET, "FIELD", order);
152  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD", order);
153  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD", order);
154 
155  }
156  if (space == HDIV) {
157  CHKERR m_field.set_field_order(root_set, MBTET, "FIELD", order);
158  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD", order);
159  }
160 
161  /****/
162  // build database
163  // build field
164  CHKERR m_field.build_fields();
165 
166  // build finite elemnts
167  CHKERR m_field.build_finite_elements();
168 
169  // build adjacencies
170  CHKERR m_field.build_adjacencies(bit_level0);
171 
172 
173  /****/
174  ProblemsManager *prb_mng_ptr;
175  CHKERR m_field.getInterface(prb_mng_ptr);
176 
177  // build problem
178  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
179 
180  // mesh partitioning
181  // partition
182  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
183  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
184 
185  // what are ghost nodes, see Petsc Manual
186  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
187 
188  typedef tee_device<std::ostream, std::ofstream> TeeDevice;
189  typedef stream<TeeDevice> TeeStream;
190 
191  std::ofstream ofs("forces_and_sources_checking_direvatives.txt");
192  TeeDevice my_tee(std::cout, ofs);
193  TeeStream my_split(my_tee);
194 
195  struct OpCheckingDirevatives
197 
198  TeeStream &mySplit;
199  OpCheckingDirevatives(TeeStream &my_split)
201  "FIELD", UserDataOperator::OPROW),
202  mySplit(my_split) {}
203 
204  MoFEMErrorCode doWork(int side, EntityType type,
207 
208  if (data.getFieldData().size() == 0)
210 
211  mySplit << "Type " << type << " side " << side << endl;
212 
213  if (data.getFieldDofs()[0]->getSpace() == H1) {
214 
215  // mySplit << std::fixed << data.getN() << std::endl;
216  // mySplit << std::fixed << data.getDiffN() << std::endl;
217 
218  const int nb_dofs = data.getN().size2();
219  for (int dd = 0; dd != nb_dofs; dd++) {
220  const double dksi = (data.getN()(1, dd) - data.getN()(0, dd)) / eps;
221  const double deta = (data.getN()(3, dd) - data.getN()(2, dd)) / eps;
222  const double dzeta =
223  (data.getN()(5, dd) - data.getN()(4, dd)) / eps;
224  mySplit << "DKsi " << dksi - data.getDiffN()(6, 3 * dd + 0)
225  << std::endl;
226  mySplit << "DEta " << deta - data.getDiffN()(6, 3 * dd + 1)
227  << std::endl;
228  mySplit << "DZeta " << dzeta - data.getDiffN()(6, 3 * dd + 2)
229  << std::endl;
230  if (fabs(dksi - data.getDiffN()(6, 3 * dd + 0)) > eps_diff) {
231  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
232  "H1 inconsistent dKsi derivative");
233  }
234  if (fabs(deta - data.getDiffN()(6, 3 * dd + 1)) > eps_diff) {
235  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
236  "H1 inconsistent dEta derivative");
237  }
238  if (fabs(dzeta - data.getDiffN()(6, 3 * dd + 2)) > eps_diff) {
239  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
240  "H1 inconsistent dZeta derivative");
241  }
242  }
243  }
244 
245  if (data.getFieldDofs()[0]->getSpace() == HDIV ||
246  data.getFieldDofs()[0]->getSpace() == HCURL) {
247 
249  &data.getN()(0, 0), &data.getN()(0, 1), &data.getN()(0, 2), 3);
251  &data.getN()(1, 0), &data.getN()(1, 1), &data.getN()(1, 2), 3);
253  &data.getN()(2, 0), &data.getN()(2, 1), &data.getN()(2, 2), 3);
255  &data.getN()(3, 0), &data.getN()(3, 1), &data.getN()(3, 2), 3);
256  FTensor::Tensor1<double *, 3> base_zeta_m(
257  &data.getN()(4, 0), &data.getN()(4, 1), &data.getN()(4, 2), 3);
258  FTensor::Tensor1<double *, 3> base_zeta_p(
259  &data.getN()(5, 0), &data.getN()(5, 1), &data.getN()(5, 2), 3);
260 
262  &data.getDiffN()(6, 0), &data.getDiffN()(6, 3),
263  &data.getDiffN()(6, 6), &data.getDiffN()(6, 1),
264  &data.getDiffN()(6, 4), &data.getDiffN()(6, 7),
265  &data.getDiffN()(6, 2), &data.getDiffN()(6, 5),
266  &data.getDiffN()(6, 8), 9);
267 
272 
273  const int nb_dofs = data.getN().size2() / 3;
274  for (int dd = 0; dd != nb_dofs; dd++) {
276  dksi(i) = (base_ksi_p(i) - base_ksi_m(i)) / eps;
278  deta(i) = (base_eta_p(i) - base_eta_m(i)) / eps;
280  dzeta(i) = (base_zeta_p(i) - base_zeta_m(i)) / eps;
281 
282  dksi(i) -= diff_base(i, N0);
283  deta(i) -= diff_base(i, N1);
284  dzeta(i) -= diff_base(i, N2);
285 
286  mySplit << "dKsi " << dksi(0) << " " << dksi(1) << " " << dksi(2)
287  << " " << sqrt(dksi(i) * dksi(i)) << endl;
288  mySplit << "dEta " << deta(0) << " " << deta(1) << " " << deta(2)
289  << " " << sqrt(deta(i) * deta(i)) << endl;
290  mySplit << "dZeta " << dzeta(0) << " " << dzeta(1) << " "
291  << dzeta(2) << " " << sqrt(dzeta(i) * dzeta(i)) << endl;
292 
293  if (sqrt(dksi(i) * dksi(i)) > eps_diff) {
294  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
295  "%s inconsistent dKsi derivative for type %d",
296  FieldSpaceNames[data.getFieldDofs()[0]->getSpace()],
297  type);
298  }
299  if (sqrt(deta(i) * deta(i)) > eps_diff) {
300  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
301  "%s inconsistent dEta derivative for type %d",
302  FieldSpaceNames[data.getFieldDofs()[0]->getSpace()],
303  type);
304  }
305  if (sqrt(dzeta(i) * dzeta(i)) > eps_diff) {
306  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
307  "%s inconsistent dZeta derivative for type %d",
308  FieldSpaceNames[data.getFieldDofs()[0]->getSpace()],
309  type);
310  }
311 
312  ++base_ksi_m;
313  ++base_ksi_p;
314  ++base_eta_m;
315  ++base_eta_p;
316  ++base_zeta_m;
317  ++base_zeta_p;
318  ++diff_base;
319  }
320  }
321 
323  }
324  };
325 
326  struct MyFE : public VolumeElementForcesAndSourcesCore {
327 
328  MyFE(MoFEM::Interface &m_field)
330  int getRule(int order) { return -1; };
331 
332  MoFEMErrorCode setGaussPts(int order) {
334 
335  const double ksi = G_TET_X1[0];
336  const double eta = G_TET_Y1[0];
337  const double zeta = G_TET_Z1[0];
338 
339  gaussPts.resize(4, 7);
340  gaussPts.clear();
341 
342  gaussPts(0, 0) = ksi - eps;
343  gaussPts(1, 0) = eta;
344  gaussPts(2, 0) = zeta;
345  gaussPts(0, 1) = ksi + eps;
346  gaussPts(1, 1) = eta;
347  gaussPts(2, 1) = zeta;
348 
349  gaussPts(0, 2) = ksi;
350  gaussPts(1, 2) = eta - eps;
351  gaussPts(2, 2) = zeta;
352  gaussPts(0, 3) = ksi;
353  gaussPts(1, 3) = eta + eps;
354  gaussPts(2, 3) = zeta;
355 
356  gaussPts(0, 4) = ksi;
357  gaussPts(1, 4) = eta;
358  gaussPts(2, 4) = zeta - eps;
359  gaussPts(0, 5) = ksi;
360  gaussPts(1, 5) = eta;
361  gaussPts(2, 5) = zeta + eps;
362 
363  gaussPts(0, 6) = ksi;
364  gaussPts(1, 6) = eta;
365  gaussPts(2, 6) = zeta;
366 
367  for (unsigned int ii = 0; ii != gaussPts.size2(); ii++) {
368  gaussPts(3, ii) = 1;
369  }
370 
371  cerr << gaussPts << endl;
372 
374  }
375  };
376 
377  MyFE tet_fe(m_field);
378 
379  tet_fe.getOpPtrVector().push_back(new OpCheckingDirevatives(my_split));
380  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TET_FE", tet_fe);
381 
382  }
383  CATCH_ERRORS;
384 
386 
387 }
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
static const double eps_diff
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
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
field with continuous normal traction
Definition: definitions.h:173
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
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
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
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.
static const double G_TET_X1[]
Definition: fem_tools.h:486
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
static const double G_TET_Y1[]
Definition: fem_tools.h:489
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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
static const double G_TET_Z1[]
Definition: fem_tools.h:492
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
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
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
int main(int argc, char *argv[])
#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
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....
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
static const double eps
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