v0.10.0
rd_stdOperators.hpp
Go to the documentation of this file.
1 #ifndef __STDOPERATORS_HPP__
2 #define __STDOPERATORS_HPP__
3 
4 #include <stdlib.h>
5 #include <BasicFiniteElements.hpp>
6 
7 namespace StdRDOperators {
8 
11 
14 
16 
17 // const double alpha = 0.01;
18 // const double gma = 0.002;
19 // const double b = 0.15;
20 // const double c = 8.00;
21 // const double mu1 = 0.20;
22 // const double mu2 = 0.30;
23 
24 const double B = 0.0;
26 const double natural_bc_values = 0.0;
27 const double essential_bc_values = 0.0;
28 const int order = 2;
29 // const int dim = 3;
31 
32 struct BlockData {
33  int block_id;
36 
37  Range block_ents;
38 
39  double B0; // species mobility
40 
41  BlockData() : B0(2e-3) {
42  coef.resize(3, 3, false);
43  rate.resize(3, false);
44  coef.clear();
45  rate.clear();
46  coef(0, 0) = 1.0;
47  coef(0, 1) = 0.0;
48  coef(0, 2) = 0.0;
49  coef(1, 0) = 0.0;
50  coef(1, 1) = 1.0;
51  coef(1, 2) = 0.0;
52  coef(2, 0) = 0.0;
53  coef(2, 1) = 0.0;
54  coef(2, 2) = 1.0;
55 
56  for (int i = 0; i < 3; ++i) {
57  rate[i] = 1.0;
58  }
59  }
60 };
61 
62 struct PreviousData {
63  MatrixDouble grads; ///< Gradients of field "u" at integration points
64  VectorDouble values; ///< Values of field "u" at integration points
66  dot_values; ///< Rate of values of field "u" at integration points
68 
69  MatrixDouble invJac; ///< Inverse of element jacobian
70 
72 };
73 
75  OpAssembleMass(std::string fieldu, SmartPetscObj<Mat> m)
76  : OpEle(fieldu, fieldu, OpEle::OPROWCOL), M(m)
77 
78  {
79  sYmm = true;
80  }
81  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
82  EntityType col_type, EntData &row_data,
83  EntData &col_data) {
85  const int nb_row_dofs = row_data.getIndices().size();
86  const int nb_col_dofs = col_data.getIndices().size();
87  if (nb_row_dofs && nb_col_dofs) {
88  const int nb_integration_pts = getGaussPts().size2();
89  mat.resize(nb_row_dofs, nb_col_dofs, false);
90  mat.clear();
91  auto t_row_base = row_data.getFTensor0N();
92  auto t_w = getFTensor0IntegrationWeight();
93  const double vol = getMeasure();
94  for (int gg = 0; gg != nb_integration_pts; ++gg) {
95  const double a = t_w * vol;
96  for (int rr = 0; rr != nb_row_dofs; ++rr) {
97  auto t_col_base = col_data.getFTensor0N(gg, 0);
98  for (int cc = 0; cc != nb_col_dofs; ++cc) {
99  mat(rr, cc) += a * t_row_base * t_col_base;
100  ++t_col_base;
101  }
102  ++t_row_base;
103  }
104  ++t_w;
105  }
106  CHKERR MatSetValues(M, row_data, col_data, &mat(0, 0), ADD_VALUES);
107  if (row_side != col_side || row_type != col_type) {
108  transMat.resize(nb_col_dofs, nb_row_dofs, false);
109  noalias(transMat) = trans(mat);
110  CHKERR MatSetValues(M, col_data, row_data, &transMat(0, 0), ADD_VALUES);
111  }
112  }
114  }
115 
116 private:
118  SmartPetscObj<Mat> M;
119 };
120 
121 struct OpComputeSlowValue : public OpEle {
122  OpComputeSlowValue(std::string mass_field,
123  boost::shared_ptr<PreviousData> &data1,
124  boost::shared_ptr<PreviousData> &data2,
125  boost::shared_ptr<PreviousData> &data3,
126  std::map<int, BlockData> &block_map)
127  : OpEle(mass_field, OpEle::OPROW), commonData1(data1), commonData2(data2),
128  commonData3(data3), massField(mass_field), setOfBlock(block_map) {}
129  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
131  boost::shared_ptr<VectorDouble> slow_value_ptr1(commonData1,
132  &commonData1->slow_values);
133  boost::shared_ptr<VectorDouble> slow_value_ptr2(commonData2,
134  &commonData2->slow_values);
135  boost::shared_ptr<VectorDouble> slow_value_ptr3(commonData3,
136  &commonData3->slow_values);
137 
138  VectorDouble &vec1 = *slow_value_ptr1;
139  VectorDouble &vec2 = *slow_value_ptr2;
140  VectorDouble &vec3 = *slow_value_ptr3;
141  const int nb_integration_pts = getGaussPts().size2();
142  if (type == MBVERTEX) {
143  vec1.resize(nb_integration_pts, false);
144  vec2.resize(nb_integration_pts, false);
145  vec3.resize(nb_integration_pts, false);
146  vec1.clear();
147  vec2.clear();
148  vec3.clear();
149  }
150  const int nb_dofs = data.getIndices().size();
151 
152  if (nb_dofs) {
153  auto find_block_data = [&]() {
154  EntityHandle fe_ent = getFEEntityHandle();
155  BlockData *block_raw_ptr = nullptr;
156  for (auto &m : setOfBlock) {
157  if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
158  block_raw_ptr = &m.second;
159  break;
160  }
161  }
162  return block_raw_ptr;
163  };
164 
165  auto block_data_ptr = find_block_data();
166  if (!block_data_ptr)
167  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
168 
169  auto &block_data = *block_data_ptr;
170 
171  const int nb_integration_pts = getGaussPts().size2();
172 
173  auto t_slow_values1 = getFTensor0FromVec(vec1);
174  auto t_slow_values2 = getFTensor0FromVec(vec2);
175  auto t_slow_values3 = getFTensor0FromVec(vec3);
176 
177  auto t_mass_values1 = getFTensor0FromVec(commonData1->values);
178  auto t_mass_values2 = getFTensor0FromVec(commonData2->values);
179  auto t_mass_values3 = getFTensor0FromVec(commonData3->values);
180 
181  // cout << "r1 : " << block_data.rate[0] << endl;
182 
183  for (int gg = 0; gg != nb_integration_pts; ++gg) {
184  t_slow_values1 = block_data.rate[0] * t_mass_values1 *
185  (1.0 - block_data.coef(0, 0) * t_mass_values1 -
186  block_data.coef(0, 1) * t_mass_values2 -
187  block_data.coef(0, 2) * t_mass_values3);
188  t_slow_values2 = block_data.rate[1] * t_mass_values2 *
189  (1.0 - block_data.coef(1, 0) * t_mass_values1 -
190  block_data.coef(1, 1) * t_mass_values2 -
191  block_data.coef(1, 2) * t_mass_values3);
192 
193  t_slow_values3 = block_data.rate[2] * t_mass_values3 *
194  (1.0 - block_data.coef(2, 0) * t_mass_values1 -
195  block_data.coef(2, 1) * t_mass_values2 -
196  block_data.coef(2, 2) * t_mass_values3);
197  ++t_slow_values1;
198  ++t_slow_values2;
199  ++t_slow_values3;
200 
201  ++t_mass_values1;
202  ++t_mass_values2;
203  ++t_mass_values3;
204  }
205  }
207  }
208 
209 private:
210  std::string massField;
211  boost::shared_ptr<PreviousData> commonData1;
212  boost::shared_ptr<PreviousData> commonData2;
213  boost::shared_ptr<PreviousData> commonData3;
214  std::map<int, BlockData> setOfBlock;
215 };
216 
218  typedef boost::function<double(const double, const double, const double)>
220  typedef boost::function<FTensor::Tensor1<double, 3>(
221  const double, const double, const double)>
223  OpAssembleSlowRhs(std::string field,
224  boost::shared_ptr<PreviousData> &data,
225  FVal exact_value,
226  FVal exact_dot,
227  FVal exact_lap)
228  : OpEle(field, OpEle::OPROW), commonData(data), exactValue(exact_value),
229  exactDot(exact_dot), exactLap(exact_lap) {}
230 
231  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
233 
234  const int nb_dofs = data.getIndices().size();
235  if (nb_dofs) {
236  vecF.resize(nb_dofs, false);
237  vecF.clear();
238  const int nb_integration_pts = getGaussPts().size2();
239 
240  auto t_slow_value = getFTensor0FromVec(commonData->slow_values);
241 
242  auto t_base = data.getFTensor0N();
243  auto t_w = getFTensor0IntegrationWeight();
244  const double vol = getMeasure();
245 
246  // cout << "measure : " << getMeasure() << endl;
247  const double ct = getFEMethod()->ts_t;
248  auto t_coords = getFTensor1CoordsAtGaussPts();
249  for (int gg = 0; gg != nb_integration_pts; ++gg) {
250  const double a = vol * t_w;
251 
252  double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
253  double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
254 
255  // double f = u_dot - u_lap;
256  const double f = t_slow_value;
257 
258  for (int rr = 0; rr != nb_dofs; ++rr) {
259  const double b = a * f * t_base;
260  vecF[rr] += b;
261  ++t_base;
262  }
263  ++t_slow_value;
264  ++t_w;
265  ++t_coords;
266  }
267  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
268  PETSC_TRUE);
269  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
270  ADD_VALUES);
271  }
273  }
274 
275 private:
276  boost::shared_ptr<PreviousData> commonData;
278 
282 
286 };
287 
288 
289 template <int DIM> struct OpAssembleStiffRhs : OpEle {
290  typedef boost::function<double(const double, const double, const double)>
292  typedef boost::function<FTensor::Tensor1<double, 3>(
293  const double, const double, const double)>
295 
296  OpAssembleStiffRhs(std::string field, boost::shared_ptr<PreviousData> &data,
297  std::map<int, BlockData> &block_map, FVal exact_value,
298  FVal exact_dot, FVal exact_lap)
299  : OpEle(field, OpEle::OPROW), commonData(data), setOfBlock(block_map),
300  exactValue(exact_value), exactDot(exact_dot), exactLap(exact_lap) {}
301 
302  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
304  const int nb_dofs = data.getIndices().size();
305  if (nb_dofs) {
306  auto find_block_data = [&]() {
307  EntityHandle fe_ent = getFEEntityHandle();
308  BlockData *block_raw_ptr = nullptr;
309  for (auto &m : setOfBlock) {
310  if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
311  block_raw_ptr = &m.second;
312  break;
313  }
314  }
315  return block_raw_ptr;
316  };
317 
318  auto block_data_ptr = find_block_data();
319  if (!block_data_ptr)
320  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
321 
322  auto &block_data = *block_data_ptr;
323 
324  vecF.resize(nb_dofs, false);
325  vecF.clear();
326  const int nb_integration_pts = getGaussPts().size2();
327 
328  auto t_dot_val = getFTensor0FromVec(commonData->dot_values);
329  auto t_val = getFTensor0FromVec(commonData->values);
330  auto t_grad = getFTensor1FromMat<DIM>(commonData->grads);
331 
332  auto t_base = data.getFTensor0N();
333  auto t_diff_base = data.getFTensor1DiffN<DIM>();
334  auto t_w = getFTensor0IntegrationWeight();
335 
337 
338  // cout << "B0 : " << block_data.B0 << endl;
339 
340  const double vol = getMeasure();
341 
342  const double ct = getFEMethod()->ts_t;
343  auto t_coords = getFTensor1CoordsAtGaussPts();
344 
345  for (int gg = 0; gg != nb_integration_pts; ++gg) {
346  const double a = vol * t_w;
347 
348  double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
349  // double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
350 
351  // double f = u_dot - u_lap;
352  double u_lap =
353  -block_data.B0 * exactLap(t_coords(NX), t_coords(NY), ct);
354  // cout << "B0 : " << block_data.B0 << endl;
355  // double f = u_dot + u_lap;
356  double f = 0.0;
357 
358  for (int rr = 0; rr != nb_dofs; ++rr) {
359  vecF[rr] +=
360  a * (t_base * t_dot_val +
361  (block_data.B0 + B * t_val) * t_diff_base(i) * t_grad(i) -
362  t_base * f);
363  ++t_diff_base;
364  ++t_base;
365  }
366  ++t_dot_val;
367  ++t_grad;
368  ++t_val;
369  ++t_w;
370  ++t_coords;
371  }
372  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
373  PETSC_TRUE);
374  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
375  ADD_VALUES);
376  }
378  }
379 
380 private:
381  boost::shared_ptr<PreviousData> commonData;
382  std::map<int, BlockData> setOfBlock;
384  std::string field;
385 
389 
393 };
394 
395 template <int DIM> struct OpAssembleStiffLhs : OpEle {
396 
397  OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr<PreviousData> &data,
398  std::map<int, BlockData> &block_map)
399  : OpEle(fieldu, fieldu, OpEle::OPROWCOL), commonData(data),
401  sYmm = true;
402  }
403  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
404  EntityType col_type, EntData &row_data,
405  EntData &col_data) {
407 
408  const int nb_row_dofs = row_data.getIndices().size();
409  const int nb_col_dofs = col_data.getIndices().size();
410  // cerr << "In doWork() : (row, col) = (" << nb_row_dofs << ", " <<
411  // nb_col_dofs << ")" << endl;
412  if (nb_row_dofs && nb_col_dofs) {
413  auto find_block_data = [&]() {
414  EntityHandle fe_ent = getFEEntityHandle();
415  BlockData *block_raw_ptr = nullptr;
416  for (auto &m : setOfBlock) {
417  if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
418  block_raw_ptr = &m.second;
419  break;
420  }
421  }
422  return block_raw_ptr;
423  };
424 
425  auto block_data_ptr = find_block_data();
426  if (!block_data_ptr)
427  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
428 
429  auto &block_data = *block_data_ptr;
430 
431  mat.resize(nb_row_dofs, nb_col_dofs, false);
432  mat.clear();
433  const int nb_integration_pts = getGaussPts().size2();
434  auto t_row_base = row_data.getFTensor0N();
435 
436  auto t_val = getFTensor0FromVec(commonData->values);
437  auto t_grad = getFTensor1FromMat<DIM>(commonData->grads);
438 
439  auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
440  auto t_w = getFTensor0IntegrationWeight();
441  const double ts_a = getFEMethod()->ts_a;
442  const double vol = getMeasure();
443 
445 
446  // cout << "B0 : " << block_data.B0 << endl;
447 
448  for (int gg = 0; gg != nb_integration_pts; ++gg) {
449  const double a = vol * t_w;
450 
451  for (int rr = 0; rr != nb_row_dofs; ++rr) {
452  auto t_col_base = col_data.getFTensor0N(gg, 0);
453  auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
454  for (int cc = 0; cc != nb_col_dofs; ++cc) {
455 
456  mat(rr, cc) +=
457  a * (t_row_base * t_col_base * ts_a +
458  (block_data.B0 + B * t_val) * t_row_diff_base(i) *
459  t_col_diff_base(i) +
460  B * t_col_base * t_grad(i) * t_row_diff_base(i));
461 
462  ++t_col_base;
463  ++t_col_diff_base;
464  }
465  ++t_row_base;
466  ++t_row_diff_base;
467  }
468  ++t_w;
469  ++t_val;
470  ++t_grad;
471  }
472  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
473  ADD_VALUES);
474  if (row_side != col_side || row_type != col_type) {
475  transMat.resize(nb_col_dofs, nb_row_dofs, false);
476  noalias(transMat) = trans(mat);
477  CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
478  &transMat(0, 0), ADD_VALUES);
479  }
480  }
482  }
483 
484  private :
485  boost::shared_ptr<PreviousData> commonData;
486  std::map<int, BlockData> setOfBlock;
488 };
489 
491 {
492  OpAssembleNaturalBCRhs(std::string mass_field, Range &natural_bd_ents)
493  : OpBoundaryEle(mass_field, OpBoundaryEle::OPROW),
495  cerr << "OpAssembleNaturalBCRhsTau()" << endl;
496  }
497 
498  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
500  const int nb_dofs = data.getIndices().size();
501 
502  if (nb_dofs) {
503  EntityHandle row_side_ent = getFEEntityHandle();
504  bool is_natural =
505  (natural_bd_ents.find(row_side_ent) != natural_bd_ents.end());
506  if (is_natural) {
507  // cerr << "In NaturalBCRhsTau..." << endl;
508  vecF.resize(nb_dofs, false);
509  vecF.clear();
510  const int nb_integration_pts = getGaussPts().size2();
511  auto t_row_base = data.getFTensor0N();
512 
513  auto dir = getDirection();
514  FTensor::Tensor1<double, 3> t_normal(-dir[1], dir[0], dir[2]);
516  auto t_w = getFTensor0IntegrationWeight();
517  const double pi = 3.141592653589793;
518  const double ct = getFEMethod()->ts_t;
519  for (int gg = 0; gg != nb_integration_pts; ++gg) {
520  const double a = t_w;
521  double x = getCoordsAtGaussPts()(gg, 0);
522  double y = getCoordsAtGaussPts()(gg, 1);
523 
524  double mm = -10 * 8 * pi * cos(2 * pi * x) * sin(2 * pi * y) *
525  sin(2 * pi * ct);
526  double nn = -10 * 8 * pi * sin(2 * pi * x) * cos(2 * pi * y) *
527  sin(2 * pi * ct);
528 
529  FTensor::Tensor1<double, 3> t_bd_val(mm, nn, 0.0);
530  double h = t_bd_val(i) * t_normal(i);
531  for (int rr = 0; rr != nb_dofs; ++rr) {
532  vecF[rr] += t_row_base * h * a;
533  ++t_row_base;
534  }
535  ++t_w;
536  }
537  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
538  PETSC_TRUE);
539  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
540  ADD_VALUES);
541  }
542  }
544  }
545 
546 private:
549 };
550 
551 struct OpError : public OpEle {
552  typedef boost::function<double(const double, const double, const double)>
554  typedef boost::function<FTensor::Tensor1<double, 3>(
555  const double, const double, const double)>
557  double &eRror;
558  OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad,
559  boost::shared_ptr<PreviousData> &prev_data,
560  std::map<int, BlockData> &block_map, double &err)
561  : OpEle("ERROR", OpEle::OPROW), exactVal(exact_value),
562  exactLap(exact_lap), exactGrad(exact_grad), prevData(prev_data),
563  setOfBlock(block_map), eRror(err) {}
564  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
566  const int nb_dofs = data.getFieldData().size();
567  // cout << "nb_error_dofs : " << nb_dofs << endl;
568  if (nb_dofs) {
569 
570  auto find_block_data = [&]() {
571  EntityHandle fe_ent = getFEEntityHandle();
572  BlockData *block_raw_ptr = nullptr;
573  for (auto &m : setOfBlock) {
574  if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
575  block_raw_ptr = &m.second;
576  break;
577  }
578  }
579  return block_raw_ptr;
580  };
581 
582  auto block_data_ptr = find_block_data();
583  if (!block_data_ptr)
584  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
585 
586  auto &block_data = *block_data_ptr;
587 
588  auto t_value = getFTensor0FromVec(prevData->values);
589  auto t_grad = getFTensor1FromMat<2>(prevData->grads);
590  // cout << "t_grad : " << t_grad << endl;
591  // auto t_grad = getFTensor1FromMat<3>(prevData->grads);
592  data.getFieldData().clear();
593  const double vol = getMeasure();
594  const int nb_integration_pts = getGaussPts().size2();
595  auto t_w = getFTensor0IntegrationWeight();
596  double dt;
597  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
598  double ct = getFEMethod()->ts_t - dt;
599  auto t_coords = getFTensor1CoordsAtGaussPts();
600 
602  // cout << "B0 : " << block_data.B0 << endl;
603 
604  for (int gg = 0; gg != nb_integration_pts; ++gg) {
605  const double a = vol * t_w;
606 
607  double mass_exact = exactVal(t_coords(NX), t_coords(NY), ct);
608  double flux_lap =
609  -block_data.B0 * exactLap(t_coords(NX), t_coords(NY), ct);
610  auto flux_exact = exactGrad(t_coords(NX), t_coords(NY), ct);
611 
612  // cout << "grad_exact : " << flux_exact << endl;
613  // cout << "grad_value : " << t_grad << endl;
614  // cout << "--------------------" << endl;
615 
616  double flux_error =
617  pow(block_data.B0, 2) * (pow(flux_exact(0) - t_grad(0), 2) +
618  pow(flux_exact(1) - t_grad(1), 2));
619 
620  double local_error = pow(mass_exact - t_value, 2); // + flux_error;
621 
622  data.getFieldData()[0] += a * local_error;
623  eRror += a * local_error;
624 
625  ++t_w;
626  ++t_value;
627  ++t_grad;
628  ++t_coords;
629  }
630 
631  data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
632  }
634  }
635 
636 private:
640  boost::shared_ptr<PreviousData> prevData;
641  std::map<int, BlockData> setOfBlock;
642 
645 };
646 
647 struct Monitor : public FEMethod {
648  double &eRror;
649  Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj<DM> &dm,
650  boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc, double &err)
651  : cOmm(comm), rAnk(rank), dM(dm), postProc(post_proc), eRror(err){};
652  MoFEMErrorCode preProcess() { return 0; }
653  MoFEMErrorCode operator()() { return 0; }
656  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-save_every_nth_step",
657  &save_every_nth_step, PETSC_NULL);
658  if (ts_step % save_every_nth_step == 0) {
660  CHKERR postProc->writeFile(
661  "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
662  }
663  Vec error_per_proc;
664  CHKERR VecCreateMPI(cOmm, 1, PETSC_DECIDE, &error_per_proc);
665  auto get_global_error = [&]() {
667  CHKERR VecSetValue(error_per_proc, rAnk, eRror, INSERT_VALUES);
669  };
670  CHKERR get_global_error();
671  CHKERR VecAssemblyBegin(error_per_proc);
672  CHKERR VecAssemblyEnd(error_per_proc);
673  double error_sum;
674  CHKERR VecSum(error_per_proc, &error_sum);
675  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Error : %3.4e \n", error_sum);
676  eRror = 0;
678  }
679 
680 private:
681  SmartPetscObj<DM> dM;
682  boost::shared_ptr<PostProcFaceOnRefinedMesh> postProc;
683  MPI_Comm cOmm;
684  const int rAnk;
685 };
686 
687 }; // namespace StdRDOperators
688 
689 #endif
static Index< 'm', 3 > m
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
#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
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
EdgeElementForcesAndSourcesCoreSwitch< 0 > EdgeElementForcesAndSourcesCore
Edge finite element default.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
double dt
Definition: heat_method.cpp:40
FTensor::Index< 'j', 3 > j
static double const pi
Definition: clipper.cpp:53
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:143
constexpr auto VecSetValues
constexpr auto MatSetValues
FTensor::Index< 'i', 3 > i
const double B
DataForcesAndSourcesCore::EntData EntData
const double essential_bc_values
FaceElementForcesAndSourcesCore::UserDataOperator OpEle
const double natural_bc_values
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
PetscReal ts_t
time
MoFEMErrorCode preProcess()
SmartPetscObj< DM > dM
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProc
MoFEMErrorCode operator()()
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcFaceOnRefinedMesh > &post_proc, double &err)
MoFEMErrorCode postProcess()
OpAssembleMass(std::string fieldu, SmartPetscObj< Mat > m)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleNaturalBCRhs(std::string mass_field, Range &natural_bd_ents)
boost::function< FTensor::Tensor1< double, 3 > const double, const double, const double)> FGrad
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleSlowRhs(std::string field, boost::shared_ptr< PreviousData > &data, FVal exact_value, FVal exact_dot, FVal exact_lap)
boost::function< double(const double, const double, const double)> FVal
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleStiffLhs(std::string fieldu, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleStiffRhs(std::string field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map, FVal exact_value, FVal exact_dot, FVal exact_lap)
boost::function< double(const double, const double, const double)> FVal
boost::function< FTensor::Tensor1< double, 3 > const double, const double, const double)> FGrad
std::map< int, BlockData > setOfBlock
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpComputeSlowValue(std::string mass_field, boost::shared_ptr< PreviousData > &data1, boost::shared_ptr< PreviousData > &data2, boost::shared_ptr< PreviousData > &data3, std::map< int, BlockData > &block_map)
std::map< int, BlockData > setOfBlock
boost::shared_ptr< PreviousData > commonData2
boost::shared_ptr< PreviousData > commonData3
boost::shared_ptr< PreviousData > commonData1
FTensor::Number< 1 > NY
boost::shared_ptr< PreviousData > prevData
std::map< int, BlockData > setOfBlock
OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad, boost::shared_ptr< PreviousData > &prev_data, std::map< int, BlockData > &block_map, double &err)
boost::function< double(const double, const double, const double)> FVal
FTensor::Number< 0 > NX
boost::function< FTensor::Tensor1< double, 3 > const double, const double, const double)> FGrad
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
VectorDouble dot_values
Rate of values of field "u" at integration points.
MatrixDouble invJac
Inverse of element jacobian.
VectorDouble values
Values of field "u" at integration points.
MatrixDouble grads
Gradients of field "u" at integration points.