v0.9.0
h1_check_approximation_2d.cpp
Go to the documentation of this file.
1 /**
2  * \file h1_check_approximation_2d.cpp
3  * \example h1_check_approximation_2d.cpp
4  *
5  * Approximate polynomial in 2D and check direvatives
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #include <MoFEM.hpp>
24 
25 using namespace MoFEM;
26 
27 static char help[] = "...\n\n";
28 
30  FaceElementForcesAndSourcesCore::NO_HO_GEOMETRY |
31  FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV |
32  FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>;
33 
36 
38  static double fUn(const double x, const double y) {
39  return pow(x, 4) + pow(y, 4) + pow(x, 2) * pow(y, 2) + x * y + x + y;
40  }
41 
42  static FTensor::Tensor1<double, 2> diffFun(const double x, const double y) {
44  4 * pow(x, 3) + 2 * x * pow(y, 2) + y + 1.,
45  4 * pow(y, 3) + 2 * y * pow(x, 2) + x + 1.);
46  }
47 };
48 
49 struct OpAssembleMat : public FaceEleOp {
50  OpAssembleMat() : FaceEleOp("FIELD1", "FIELD1", OPROWCOL) {
51  sYmm = false; // FIXME problem is symmetric, should use that
52  }
53 
55  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
56  EntityType col_type, EntData &row_data,
57  EntData &col_data) {
58 
60  const int nb_dofs_row = row_data.getIndices().size();
61  if (nb_dofs_row == 0)
63  const int nb_dofs_col = col_data.getIndices().size();
64  if (nb_dofs_col == 0)
66  nA.resize(nb_dofs_row, nb_dofs_col, false);
67  nA.clear();
68  const int nb_gauss_pts = row_data.getN().size1();
69  auto t_base_row = row_data.getFTensor0N();
70  for (int gg = 0; gg != nb_gauss_pts; gg++) {
71  const double val = getArea() * getGaussPts()(2, gg);
72  for (int rr = 0; rr != nb_dofs_row; rr++) {
73  auto t_base_col = col_data.getFTensor0N(gg, 0);
74  for (int cc = 0; cc != nb_dofs_col; cc++) {
75  nA(rr, cc) += val * t_base_row * t_base_col;
76  ++t_base_col;
77  }
78  ++t_base_row;
79  }
80  }
81  CHKERR MatSetValues(getFEMethod()->ksp_B, row_data, col_data,
82  &*nA.data().begin(), ADD_VALUES);
84  }
85 };
86 
87 struct OpAssembleVec : public FaceEleOp {
88  OpAssembleVec() : FaceEleOp("FIELD1", "FIELD1", OPROW) {}
89 
91 
93  MoFEMErrorCode doWork(int side, EntityType type,
95 
97  const int nb_dofs = data.getIndices().size();
98  if (nb_dofs == 0)
100  const int nb_gauss_pts = data.getN().size1();
101  nF.resize(nb_dofs, false);
102  nF.clear();
103  auto t_base_fun = data.getFTensor0N();
104  for (int gg = 0; gg != nb_gauss_pts; gg++) {
105  const double val = getArea() * getGaussPts()(2, gg);
106  for (int bb = 0; bb != nb_dofs; bb++) {
107  const double x = getCoordsAtGaussPts()(gg, 0);
108  const double y = getCoordsAtGaussPts()(gg, 1);
109  nF(bb) += val * t_base_fun * ApproxFunctions::fUn(x, y);
110  ++t_base_fun;
111  }
112  }
113  CHKERR VecSetValues(getFEMethod()->ksp_f, data, &*nF.data().begin(),
114  ADD_VALUES);
116  }
117 };
118 
119 struct OpValsDiffVals : public FaceEleOp {
123  : FaceEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals) {}
124 
126 
127  MoFEMErrorCode doWork(int side, EntityType type,
130  const int nb_dofs = data.getIndices().size();
131  if (nb_dofs == 0)
133  const int nb_gauss_pts = data.getN().size1();
134  if (type == MBVERTEX) {
135  vAls.resize(nb_gauss_pts, false);
136  diffVals.resize(2, nb_gauss_pts, false);
137  vAls.clear();
138  diffVals.clear();
139  }
140  auto t_vals = getFTensor0FromVec(vAls);
141  auto t_diff_vals = getFTensor1FromMat<2>(diffVals);
142  auto t_base_fun = data.getFTensor0N();
143  auto t_diff_base_fun = data.getFTensor1DiffN<2>();
144  for (int gg = 0; gg != nb_gauss_pts; gg++) {
145  auto t_data = data.getFTensor0FieldData();
146  for (int bb = 0; bb != nb_dofs; bb++) {
147  t_vals += t_base_fun * t_data;
148  t_diff_vals(i) += t_diff_base_fun(i) * t_data;
149  ++t_base_fun;
150  ++t_diff_base_fun;
151  ++t_data;
152  }
153  ++t_vals;
154  ++t_diff_vals;
155  }
157  }
158 };
159 
163  boost::shared_ptr<VectorDouble> ptrVals;
164  boost::shared_ptr<MatrixDouble> ptrDiffVals;
165 
167  boost::shared_ptr<VectorDouble> &ptr_vals,
168  boost::shared_ptr<MatrixDouble> &ptr_diff_vals)
169  : FaceEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
170  ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals) {}
171 
173 
174  MoFEMErrorCode doWork(int side, EntityType type,
177  const double eps = 1e-6;
178  if (type == MBVERTEX) {
179  const int nb_gauss_pts = data.getN().size1();
180 
181  auto t_vals = getFTensor0FromVec(vAls);
182  auto t_diff_vals = getFTensor1FromMat<2>(diffVals);
183 
184  auto t_ptr_vals = getFTensor0FromVec(*ptrVals);
185  auto t_ptr_diff_vals = getFTensor1FromMat<2>(*ptrDiffVals);
186 
187  for (int gg = 0; gg != nb_gauss_pts; gg++) {
188  const double x = getCoordsAtGaussPts()(gg, 0);
189  const double y = getCoordsAtGaussPts()(gg, 1);
190 
191  // Check approximation
192  const double delta_val = t_vals - ApproxFunctions::fUn(x, y);
193  FTensor::Tensor1<double, 2> t_delta_diff_val;
194  t_delta_diff_val(i) =
195  t_diff_vals(i) - ApproxFunctions::diffFun(x, y)(i);
196  double err_val = sqrt(delta_val * delta_val);
197  double err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
198  if (err_val > eps)
199  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
200  "Wrong value %4.3e", err_val);
201 
202  if (err_diff_val > eps)
203  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
204  "Wrong derivative of value %4.3e", err_diff_val);
205 
206  // Check H1 user data operators
207  err_val = t_vals - t_ptr_vals;
208  if (err_val > eps)
209  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
210  "Wrong value from operator %4.3e", err_val);
211 
212  t_delta_diff_val(i) = t_diff_vals(i) - t_ptr_diff_vals(i);
213  err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
214  if (err_diff_val > eps)
215  SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
216  "Wrong direvatives from operator %4.3e", err_diff_val);
217 
218  ++t_vals;
219  ++t_diff_vals;
220  ++t_ptr_vals;
221  ++t_ptr_diff_vals;
222  }
223  }
225  }
226 };
227 
228 int main(int argc, char *argv[]) {
229 
230  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
231 
232  try {
233 
234  DMType dm_name = "DMMOFEM";
235  CHKERR DMRegister_MoFEM(dm_name);
236 
237  moab::Core mb_instance;
238  moab::Interface &moab = mb_instance;
239 
240  // Create MoFEM instance
241  MoFEM::Core core(moab);
242  MoFEM::Interface &m_field = core;
243 
244  Simple *simple_interface = m_field.getInterface<Simple>();
245  Basic *basic_interface = m_field.getInterface<Basic>();
246  CHKERR simple_interface->getOptions();
247  CHKERR simple_interface->loadFile("", "rectangle.h5m");
248 
249  // Declare elements
250  enum bases { AINSWORTH, DEMKOWICZ, BERNSTEIN, LASBASETOP };
251  const char *list_bases[] = {"ainsworth", "demkowicz", "bernstein"};
252  PetscBool flg;
253  PetscInt choice_base_value = AINSWORTH;
254  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
255  LASBASETOP, &choice_base_value, &flg);
256 
257  if (flg != PETSC_TRUE)
258  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
260  if (choice_base_value == AINSWORTH)
262  else if (choice_base_value == DEMKOWICZ)
263  base = DEMKOWICZ_JACOBI_BASE;
264  else if (choice_base_value == BERNSTEIN)
266 
267  CHKERR simple_interface->addDomainField("FIELD1", H1, base, 1);
268  constexpr int order = 5;
269  CHKERR simple_interface->setFieldOrder("FIELD1", order);
270  CHKERR simple_interface->setUp();
271  auto dm = simple_interface->getDM();
272 
273  VectorDouble vals;
274  MatrixDouble jac(2, 2), inv_jac(2, 2), diff_vals;
275 
276  auto assemble_matrices_and_vectors = [&]() {
278  basic_interface->getOpDomainRhsPipeline().push_back(new OpAssembleVec());
279 
280  basic_interface->getOpDomainLhsPipeline().push_back(
281  new OpCalculateInvJacForFace(inv_jac));
282  basic_interface->getOpDomainLhsPipeline().push_back(new OpAssembleMat());
283 
284  auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
285  CHKERR basic_interface->setDomainRhsIntegrationRule(integration_rule);
286  CHKERR basic_interface->setDomainLhsIntegrationRule(integration_rule);
287 
289  };
290 
291  auto solve_problem = [&] {
293  auto solver = basic_interface->createKSP();
294  CHKERR KSPSetFromOptions(solver);
295  CHKERR KSPSetUp(solver);
296 
297  auto D = smartCreateDMDVector(dm);
298  auto F = smartVectorDuplicate(D);
299 
300  CHKERR KSPSolve(solver, F, D);
301  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
302  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
303  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
305  };
306 
307  auto check_solution = [&] {
309 
310  boost::shared_ptr<VectorDouble> ptr_values(new VectorDouble());
311  boost::shared_ptr<MatrixDouble> ptr_diff_vals(new MatrixDouble());
312 
313  basic_interface->getDomainLhsFE().reset();
314  basic_interface->getOpDomainRhsPipeline().clear();
315 
316  basic_interface->getOpDomainRhsPipeline().push_back(
317  new OpCalculateInvJacForFace(inv_jac));
318  basic_interface->getOpDomainRhsPipeline().push_back(
319  new OpSetInvJacH1ForFace(inv_jac));
320  basic_interface->getOpDomainRhsPipeline().push_back(
321  new OpValsDiffVals(vals, diff_vals));
322  basic_interface->getOpDomainRhsPipeline().push_back(
323  new OpCalculateScalarFieldValues("FIELD1", ptr_values));
324  basic_interface->getOpDomainRhsPipeline().push_back(
325  new OpCalculateScalarFieldGradient<2>("FIELD1", ptr_diff_vals));
326  basic_interface->getOpDomainRhsPipeline().push_back(
327  new OpCheckValsDiffVals(vals, diff_vals, ptr_values, ptr_diff_vals));
328 
329  CHKERR basic_interface->loopFiniteElements();
330 
332  };
333 
334  CHKERR assemble_matrices_and_vectors();
336  CHKERR check_solution();
337  }
338  CATCH_ERRORS;
339 
341 }
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:433
static FTensor::Tensor1< double, 2 > diffFun(const double x, const double y)
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:443
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:434
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
static char help[]
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Index< 'i', 2 > i
Core (interface) class.
Definition: Core.hpp:50
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: AuxPETSc.hpp:238
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:512
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
int main(int argc, char *argv[])
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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
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)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Basic interface.
Definition: Basic.hpp:36
Calculate inverse of jacobian for face element.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
FieldApproximationBase
approximation base
Definition: definitions.h:148
MoFEMErrorCode loadFile(const std::string options="PARALLEL=READ_PART;" "PARALLEL_RESOLVE_SHARED_ENTS;" "PARTITION=PARALLEL_PARTITION;", const std::string mesh_file_name="")
Load mesh file.
Definition: Simple.cpp:61
static double fUn(const double x, const double y)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:53
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:150
MoFEMErrorCode addDomainField(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_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:105
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:600
Transform local reference derivatives of shape functions to global derivatives.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
DataForcesAndSourcesCore::EntData EntData
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:47
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:437
continuous field
Definition: definitions.h:175
boost::shared_ptr< VectorDouble > ptrVals
static const double eps
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
auto smartCreateDMDVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:930
FTensor::Index< 'i', 2 > i
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:249
Get value at integration points for scalar field.
boost::shared_ptr< MatrixDouble > ptrDiffVals
OpCheckValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, boost::shared_ptr< VectorDouble > &ptr_vals, boost::shared_ptr< MatrixDouble > &ptr_diff_vals)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61