v0.13.0
PoissonOperators.hpp
Go to the documentation of this file.
1 /**
2  * \file PoissonOperators.hpp
3  * \example PoissonOperators.hpp
4  *
5  */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifndef __POISSONOPERATORS_HPP__
22 #define __POISSONOPERATORS_HPP__
23 
24 namespace PoissonExample {
25 
26 /**
27  * \brief Calculate the grad-grad operator and assemble matrix
28  *
29  * Calculate
30  * \f[
31  * \mathbf{K}=\int_\Omega \nabla \boldsymbol\phi \cdot \nabla \boldsymbol\phi
32  * \textrm{d}\Omega \f] and assemble to global matrix.
33  *
34  * This operator is executed on element for each unique combination of entities.
35  *
36  */
38 
39  OpK(bool symm = true)
41  symm) {}
42 
43  /**
44  * \brief Do calculations for give operator
45  * @param row_side row side number (local number) of entity on element
46  * @param col_side column side number (local number) of entity on element
47  * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
48  * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
49  * @param row_data data for row
50  * @param col_data data for column
51  * @return error code
52  */
53  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
54  EntityType col_type,
56  EntitiesFieldData::EntData &col_data) {
58  // get number of dofs on row
59  nbRows = row_data.getIndices().size();
60  // if no dofs on row, exit that work, nothing to do here
61  if (!nbRows)
63  // get number of dofs on column
64  nbCols = col_data.getIndices().size();
65  // if no dofs on Columbia, exit nothing to do here
66  if (!nbCols)
68  // get number of integration points
69  nbIntegrationPts = getGaussPts().size2();
70  // check if entity block is on matrix diagonal
71  if (row_side == col_side && row_type == col_type) {
72  isDiag = true; // yes, it is
73  } else {
74  isDiag = false;
75  }
76  // integrate local matrix for entity block
77  CHKERR iNtegrate(row_data, col_data);
78  // assemble local matrix
79  CHKERR aSsemble(row_data, col_data);
81  }
82 
83 protected:
84  ///< error code
85 
86  int nbRows; ///< number of dofs on rows
87  int nbCols; ///< number if dof on column
88  int nbIntegrationPts; ///< number of integration points
89  bool isDiag; ///< true if this block is on diagonal
90 
91  FTensor::Index<'i', 3> i; ///< summit Index
92  MatrixDouble locMat; ///< local entity block matrix
93 
94  /**
95  * \brief Integrate grad-grad operator
96  * @param row_data row data (consist base functions on row entity)
97  * @param col_data column data (consist base functions on column entity)
98  * @return error code
99  */
100  virtual MoFEMErrorCode
102  EntitiesFieldData::EntData &col_data) {
104  // set size of local entity bock
105  locMat.resize(nbRows, nbCols, false);
106  // clear matrix
107  locMat.clear();
108  // get element volume
109  double vol = getVolume();
110  // get integration weights
111  auto t_w = getFTensor0IntegrationWeight();
112  // get base function gradient on rows
113  auto t_row_grad = row_data.getFTensor1DiffN<3>();
114  // loop over integration points
115  for (int gg = 0; gg != nbIntegrationPts; gg++) {
116  // take into account Jacobean
117  const double alpha = t_w * vol;
118  // noalias(locMat) +=
119  // alpha*prod(row_data.getDiffN(gg),trans(col_data.getDiffN(gg))); take
120  // fist element to local matrix
121  FTensor::Tensor0<double *> a(&*locMat.data().begin());
122  // loop over rows base functions
123  for (int rr = 0; rr != nbRows; rr++) {
124  // get column base functions gradient at gauss point gg
125  auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
126  // loop over columns
127  for (int cc = 0; cc != nbCols; cc++) {
128  // calculate element of local matrix
129  a += alpha * (t_row_grad(i) * t_col_grad(i));
130  ++t_col_grad; // move to another gradient of base function on column
131  ++a; // move to another element of local matrix in column
132  }
133  ++t_row_grad; // move to another element of gradient of base function on
134  // row
135  }
136  ++t_w; // move to another integration weight
137  }
139  }
140 
141  /**
142  * \brief Assemble local entity block matrix
143  * @param row_data row data (consist base functions on row entity)
144  * @param col_data column data (consist base functions on column entity)
145  * @return error code
146  */
148  EntitiesFieldData::EntData &col_data) {
150  // get pointer to first global index on row
151  const int *row_indices = &*row_data.getIndices().data().begin();
152  // get pointer to first global index on column
153  const int *col_indices = &*col_data.getIndices().data().begin();
154  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
155  : getFEMethod()->snes_B;
156  // assemble local matrix
157  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
158  &*locMat.data().begin(), ADD_VALUES);
159 
160  if (!isDiag && sYmm) {
161  // if not diagonal term and since global matrix is symmetric assemble
162  // transpose term.
163  locMat = trans(locMat);
164  CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
165  &*locMat.data().begin(), ADD_VALUES);
166  }
168  }
169 };
170 
171 /**
172  * \brief template class for integration oh the right hand side
173  */
174 template <typename OPBASE> struct OpBaseRhs : public OPBASE {
175 
176  OpBaseRhs(const std::string field_name) : OPBASE(field_name, OPBASE::OPROW) {}
177 
178  /**
179  * \brief This function is called by finite element
180  *
181  * Do work is composed from two operations, integrate and assembly. Also,
182  * it set values nbRows, and nbIntegrationPts.
183  *
184  */
185  MoFEMErrorCode doWork(int row_side, EntityType row_type,
186  EntitiesFieldData::EntData &row_data) {
188  // get number of dofs on row
189  nbRows = row_data.getIndices().size();
190  if (!nbRows)
192  // get number of integration points
193  nbIntegrationPts = OPBASE::getGaussPts().size2();
194  // integrate local vector
195  CHKERR iNtegrate(row_data);
196  // assemble local vector
197  CHKERR aSsemble(row_data);
199  }
200 
201  /**
202  * \brief Class dedicated to integrate operator
203  * @param data entity data on element row
204  * @return error code
205  */
207 
208  /**
209  * \brief Class dedicated to assemble operator to global system vector
210  * @param data entity data (indices, base functions, etc. ) on element row
211  * @return error code
212  */
214 
215 protected:
216  ///< error code
217  int nbRows; ///< number of dofs on row
218  int nbIntegrationPts; ///< number of integration points
219 };
220 
221 /**
222  * \brief Operator calculate source term,
223  *
224  * \f[
225  * \mathbf{F} = \int_\Omega \boldsymbol\phi f \textrm{d}\Omega
226  * \f]
227  *
228  */
229 struct OpF
230  : public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
231 
232  typedef boost::function<double(const double, const double, const double)>
234 
235  OpF(FSource f_source)
237  fSource(f_source) {}
238 
239 protected:
244 
246 
247  /**
248  * \brief Integrate local entity vector
249  * @param data entity data on element row
250  * @return error code
251  */
254  // set size of local vector
255  locVec.resize(nbRows, false);
256  // clear local entity vector
257  locVec.clear();
258  // get finite element volume
259  double vol = getVolume();
260  // get integration weights
261  auto t_w = getFTensor0IntegrationWeight();
262  // get base functions on entity
263  auto t_v = data.getFTensor0N();
264  // get coordinates at integration points
265  auto t_coords = getFTensor1CoordsAtGaussPts();
266  // loop over all integration points
267  for (int gg = 0; gg != nbIntegrationPts; gg++) {
268  // evaluate constant term
269  const double alpha =
270  vol * t_w * fSource(t_coords(NX), t_coords(NY), t_coords(NZ));
271  // get element of local vector
273  &*locVec.data().begin());
274  // loop over base functions
275  for (int rr = 0; rr != nbRows; rr++) {
276  // add to local vector source term
277  t_a -= alpha * t_v;
278  ++t_a; // move to next element of local vector
279  ++t_v; // move to next base function
280  }
281  ++t_w; // move to next integration weight
282  ++t_coords; // move to next physical coordinates at integration point
283  }
285  }
286 
287  /**
288  * \brief assemble local entity vector to the global right hand side
289  * @param data entity data, i.e. global indices of local vector
290  * @return error code
291  */
294  // get global indices of local vector
295  const int *indices = &*data.getIndices().data().begin();
296  // get values from local vector
297  const double *vals = &*locVec.data().begin();
298  Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
299  : getFEMethod()->snes_f;
300  // assemble vector
301  CHKERR VecSetValues(f, nbRows, indices, vals, ADD_VALUES);
303  }
304 };
305 
306 /**
307  * \brief Calculate constrains matrix
308  *
309  * \f[
310  * \mathbf{C} = \int_{\partial\Omega} \boldsymbol\psi \boldsymbol\phi
311  * \textrm{d}\partial\Omega \f] where \f$\lambda \f$ is base function on
312  * boundary
313  *
314  */
316 
317  OpC(const bool assemble_transpose)
319  false),
320  assembleTranspose(assemble_transpose) {}
321 
322  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
323  EntityType col_type,
324  EntitiesFieldData::EntData &row_data,
325  EntitiesFieldData::EntData &col_data) {
327  // get number of dofs on row
328  nbRows = row_data.getIndices().size();
329  // exit here if no dofs on row, nothing to do
330  if (!nbRows)
332  // get number of dofs on column,
333  nbCols = col_data.getIndices().size();
334  // exit here if no dofs on roe, nothing to do
335  if (!nbCols)
337  // get number of integration points
338  nbIntegrationPts = getGaussPts().size2();
339  // integrate local constrains matrix
340  CHKERR iNtegrate(row_data, col_data);
341  // assemble local constrains matrix
342  CHKERR aSsemble(row_data, col_data);
344  }
345 
346 private:
347  ///< error code
348 
349  int nbRows; ///< number of dofs on row
350  int nbCols; ///< number of dofs on column
351  int nbIntegrationPts; ///< number of integration points
352  const bool assembleTranspose; ///< assemble transpose, i.e. CT if set to true
353 
354  MatrixDouble locMat; ///< local constrains matrix
355 
356  /** \brief Integrate local constrains matrix
357  */
359  EntitiesFieldData::EntData &col_data) {
361  // set size of local constrains matrix
362  locMat.resize(nbRows, nbCols, false);
363  // clear matrix
364  locMat.clear();
365  // get area of element
366  const double area = getArea();
367  // get integration weights
368  auto t_w = getFTensor0IntegrationWeight();
369  // get base functions on entity
370  auto t_row = row_data.getFTensor0N();
371  // run over integration points
372  for (int gg = 0; gg != nbIntegrationPts; gg++) {
373  const double alpha = area * t_w;
374  // get element of local matrix
376  &*locMat.data().begin());
377  // run over base functions on rows
378  for (int rr = 0; rr != nbRows; rr++) {
379  // get first base functions on column for integration point gg
380  auto t_col = col_data.getFTensor0N(gg, 0);
381  // run over base function on column
382  for (int cc = 0; cc != nbCols; cc++) {
383  // integrate element of constrains matrix
384  c += alpha * t_row * t_col;
385  ++t_col; // move to next base function on column
386  ++c; // move to next element of local matrix
387  }
388  ++t_row; // move to next base function on row
389  }
390  ++t_w; // move to next integrate weight
391  }
393  }
394 
395  /**
396  * \brief integrate local constrains matrix
397  */
399  EntitiesFieldData::EntData &col_data) {
401  // get indices on row
402  const int *row_indices = &*row_data.getIndices().data().begin();
403  // get indices on column
404  const int *col_indices = &*col_data.getIndices().data().begin();
405  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
406  : getFEMethod()->snes_B;
407  // assemble local matrix
408  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
409  &*locMat.data().begin(), ADD_VALUES);
410  // cerr << locMat << endl;
411  if (assembleTranspose) {
412  // assemble transpose of local matrix
413  locMat = trans(locMat);
414  CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
415  &*locMat.data().begin(), ADD_VALUES);
416  }
418  }
419 };
420 
421 /**
422  * \brief Assemble constrains vector
423  *
424  * \f[
425  * \mathbf{g} = \int_{\partial\Omega} \boldsymbol\psi \overline{u}
426  * \textrm{d}\partial\Omega \f]
427  *
428  */
429 struct Op_g
430  : public OpBaseRhs<FaceElementForcesAndSourcesCore::UserDataOperator> {
431 
432  typedef boost::function<double(const double, const double, const double)>
434 
435  Op_g(FVal f_value, const string field_name = "L", const double beta = 1)
437  field_name),
438  fValue(f_value), bEta(beta) {}
439 
440 protected:
441  FTensor::Number<0> NX; ///< x-direction index
442  FTensor::Number<1> NY; ///< y-direction index
443  FTensor::Number<2> NZ; ///< z-direction index
444  FVal fValue; ///< Function pointer evaluating values of "U" at the boundary
445 
447  const double bEta;
448 
449  /**
450  * \brief Integrate local constrains vector
451  */
454  // set size to local vector
455  locVec.resize(nbRows, false);
456  // clear local vector
457  locVec.clear();
458  // get face area
459  const double area = getArea() * bEta;
460  // get integration weight
461  auto t_w = getFTensor0IntegrationWeight();
462  // get base function
463  auto t_l = data.getFTensor0N();
464  // get coordinates at integration point
465  auto t_coords = getFTensor1CoordsAtGaussPts();
466  // make loop over integration points
467  for (int gg = 0; gg != nbIntegrationPts; gg++) {
468  // evaluate function on boundary and scale it by area and integration
469  // weight
470  double alpha =
471  area * t_w * fValue(t_coords(NX), t_coords(NY), t_coords(NZ));
472  // get element of vector
474  &*locVec.data().begin());
475  //
476  for (int rr = 0; rr != nbRows; rr++) {
477  t_a += alpha * t_l;
478  ++t_a;
479  ++t_l;
480  }
481  ++t_w;
482  ++t_coords;
483  }
485  }
486 
487  /**
488  * \brief assemble constrains vectors
489  */
492  const int *indices = &*data.getIndices().data().begin();
493  const double *vals = &*locVec.data().begin();
494  Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
495  : getFEMethod()->snes_f;
496  CHKERR VecSetValues(f, nbRows, indices, &*vals, ADD_VALUES);
498  }
499 };
500 
501 /**
502  * \brief Evaluate error
503  */
504 struct OpError
505  : public OpBaseRhs<VolumeElementForcesAndSourcesCore::UserDataOperator> {
506 
507  typedef boost::function<double(const double, const double, const double)>
509  typedef boost::function<FTensor::Tensor1<double, 3>(
510  const double, const double, const double)>
512 
513  OpError(UVal u_value, GVal g_value,
514  boost::shared_ptr<VectorDouble> &field_vals,
515  boost::shared_ptr<MatrixDouble> &grad_vals, Vec global_error)
517  globalError(global_error), uValue(u_value), gValue(g_value),
518  fieldVals(field_vals), gradVals(grad_vals) {}
519 
520  MoFEMErrorCode doWork(int row_side, EntityType row_type,
521  EntitiesFieldData::EntData &row_data) {
523  nbRows = row_data.getFieldData().size();
524  if (!nbRows)
526  nbIntegrationPts = getGaussPts().size2();
527  CHKERR iNtegrate(row_data);
528  CHKERR aSsemble(row_data);
530  }
531 
532 private:
533  Vec globalError; ///< ghost vector with global (integrated over volume) error
534 
539  UVal uValue; ///< function with exact solution
540  GVal gValue; ///< function with exact solution for gradient
541 
542  boost::shared_ptr<VectorDouble> fieldVals;
543  boost::shared_ptr<MatrixDouble> gradVals;
544 
545  /**
546  * \brief Integrate error
547  */
550  // clear field dofs
551  data.getFieldData().clear();
552  // get volume of element
553  const double vol = getVolume();
554  // get integration weight
555  auto t_w = getFTensor0IntegrationWeight();
556  // get solution at integration point
557  auto t_u = getFTensor0FromVec(*fieldVals);
558  // get solution at integration point
559  auto t_grad = getFTensor1FromMat<3>(*gradVals);
560  // get coordinates at integration point
561  auto t_coords = getFTensor1CoordsAtGaussPts();
562  // keep exact gradient and error or gradient
563  FTensor::Tensor1<double, 3> t_exact_grad, t_error_grad;
564  // integrate over
565  for (int gg = 0; gg != nbIntegrationPts; gg++) {
566  double alpha = vol * t_w;
567  // evaluate exact value
568  double exact_u = uValue(t_coords(NX), t_coords(NY), t_coords(NZ));
569  // evaluate exact gradient
570  t_exact_grad = gValue(t_coords(NX), t_coords(NY), t_coords(NZ));
571  // calculate gradient errro
572  t_error_grad(i) = t_grad(i) - t_exact_grad(i);
573  // error
574  double error = pow(t_u - exact_u, 2) + t_error_grad(i) * t_error_grad(i);
575  // iterate over base functions
576  data.getFieldData()[0] += alpha * error;
577  ++t_w; // move to next integration point
578  ++t_u; // next value of function at integration point
579  ++t_grad; // next gradient at integration point
580  ++t_coords; // next coordinate at integration point
581  }
583  }
584 
585  /**
586  * \brief Assemble error
587  */
590  // set error on mesh
591  data.getFieldDofs()[0]->getFieldData() = sqrt(data.getFieldData()[0]);
592  // assemble vector to global error
593  CHKERR VecSetValue(globalError, 0, data.getFieldData()[0], ADD_VALUES);
595  }
596 };
597 
598 struct OpKt : public OpK {
599 
600  OpKt(boost::function<double(const double)> a,
601  boost::function<double(const double)> diff_a,
602  boost::shared_ptr<VectorDouble> &field_vals,
603  boost::shared_ptr<MatrixDouble> &grad_vals)
604  : OpK(false), A(a), diffA(diff_a), fieldVals(field_vals),
605  gradVals(grad_vals) {}
606 
607 protected:
608  /**
609  * \brief Integrate grad-grad operator
610  * @param row_data row data (consist base functions on row entity)
611  * @param col_data column data (consist base functions on column entity)
612  * @return error code
613  */
615  EntitiesFieldData::EntData &col_data) {
617  // set size of local entity bock
618  locMat.resize(nbRows, nbCols, false);
619  // clear matrix
620  locMat.clear();
621  // get element volume
622  double vol = getVolume();
623  // get integration weights
624  auto t_w = getFTensor0IntegrationWeight();
625  // get solution at integration point
626  auto t_u = getFTensor0FromVec(*fieldVals);
627  // get solution at integration point
628  auto t_grad = getFTensor1FromMat<3>(*gradVals);
629  // get base function gradient on rows
630  auto t_row_grad = row_data.getFTensor1DiffN<3>();
631  // loop over integration points
632  for (int gg = 0; gg != nbIntegrationPts; gg++) {
633  // take into account Jacobian
634  const double alpha = t_w * vol;
635  const double beta = alpha * A(t_u);
637  t_gamma(i) = (alpha * diffA(t_u)) * t_grad(i);
638  // take fist element to local matrix
640  &*locMat.data().begin());
641  // loop over rows base functions
642  for (int rr = 0; rr != nbRows; rr++) {
643  // get column base function
644  auto t_col = col_data.getFTensor0N(gg, 0);
645  // get column base functions gradient at gauss point gg
646  auto t_col_grad = col_data.getFTensor1DiffN<3>(gg, 0);
647  // loop over columns
648  for (int cc = 0; cc != nbCols; cc++) {
649  // calculate element of local matrix
650  a += (t_row_grad(i) * beta) * t_col_grad(i) +
651  t_row_grad(i) * (t_gamma(i) * t_col);
652  ++t_col; // move to next base function
653  ++t_col_grad; // move to another gradient of base function on column
654  ++a; // move to another element of local matrix in column
655  }
656  ++t_row_grad; // move to another element of gradient of base function on
657  // row
658  }
659  ++t_w; // move to another integration weight
660  ++t_u; // move to next value at integration point
661  ++t_grad; // move to next gradient value
662  }
664  }
665 
666  boost::function<double(const double)> A;
667  boost::function<double(const double)> diffA;
668  boost::shared_ptr<VectorDouble> fieldVals;
669  boost::shared_ptr<MatrixDouble> gradVals;
670 };
671 
672 struct OpResF_Domain : public OpF {
673 
674  OpResF_Domain(FSource f_source, boost::function<double(const double)> a,
675  boost::shared_ptr<VectorDouble> &field_vals,
676  boost::shared_ptr<MatrixDouble> &grad_vals)
677  : OpF(f_source), A(a), fieldVals(field_vals), gradVals(grad_vals) {}
678 
679 protected:
681 
682  /**
683  * \brief Integrate local entity vector
684  * @param data entity data on element row
685  * @return error code
686  */
689  // set size of local vector
690  locVec.resize(nbRows, false);
691  // clear local entity vector
692  locVec.clear();
693  // get finite element volume
694  double vol = getVolume();
695  // get integration weights
696  auto t_w = getFTensor0IntegrationWeight();
697  // get solution at integration point
698  auto t_u = getFTensor0FromVec(*fieldVals);
699  // get solution at integration point
700  auto t_grad = getFTensor1FromMat<3>(*gradVals);
701  // get base functions on entity
702  auto t_v = data.getFTensor0N();
703  // get base function gradient on rows
704  auto t_v_grad = data.getFTensor1DiffN<3>();
705  // get coordinates at integration points
706  auto t_coords = getFTensor1CoordsAtGaussPts();
707  // loop over all integration points
708  for (int gg = 0; gg != nbIntegrationPts; gg++) {
709  // evaluate constant term
710  const double alpha = vol * t_w;
711  const double source_term =
712  alpha * fSource(t_coords(NX), t_coords(NY), t_coords(NZ));
713  FTensor::Tensor1<double, 3> grad_term;
714  grad_term(i) = (alpha * A(t_u)) * t_grad(i);
715  // get element of local vector
717  &*locVec.data().begin());
718  // loop over base functions
719  for (int rr = 0; rr != nbRows; rr++) {
720  // add to local vector source term
721  t_a += t_v_grad(i) * grad_term(i) + t_v * source_term;
722  ++t_a; // move to next element of local vector
723  ++t_v; // move to next base function
724  ++t_v_grad; // move to next gradient of base function
725  }
726  ++t_w; // move to next integration weights
727  ++t_u; // move to next value
728  ++t_grad; // move to next gradient value
729  ++t_coords; // move to next physical coordinates at integration point
730  }
732  }
733 
734  boost::function<double(const double)> A;
735  boost::shared_ptr<VectorDouble> fieldVals;
736  boost::shared_ptr<MatrixDouble> gradVals;
737 };
738 
739 struct OpRes_g : public Op_g {
740 
741  OpRes_g(FVal f_value, boost::shared_ptr<VectorDouble> &field_vals)
742  : Op_g(f_value, "L", 1), fieldVals(field_vals) {}
743 
744 protected:
745  boost::shared_ptr<VectorDouble> fieldVals;
746 
747  /**
748  * \brief Integrate local constrains vector
749  */
752  // set size to local vector
753  locVec.resize(nbRows, false);
754  // clear local vector
755  locVec.clear();
756  // get face area
757  const double area = getArea() * bEta;
758  // get integration weight
759  auto t_w = getFTensor0IntegrationWeight();
760  // get base function
761  auto t_l = data.getFTensor0N();
762  // get solution at integration point
763  auto t_u = getFTensor0FromVec(*fieldVals);
764  // get coordinates at integration point
765  auto t_coords = getFTensor1CoordsAtGaussPts();
766  // make loop over integration points
767  for (int gg = 0; gg != nbIntegrationPts; gg++) {
768  // evaluate function on boundary and scale it by area and integration
769  // weight
770  double alpha = area * t_w;
771  // get element of vector
773  &*locVec.data().begin());
774  for (int rr = 0; rr != nbRows; rr++) {
775  t_a += alpha * t_l *
776  (t_u - fValue(t_coords(NX), t_coords(NY), t_coords(NZ)));
777  ++t_a;
778  ++t_l;
779  }
780  ++t_w;
781  ++t_u;
782  ++t_coords;
783  }
785  }
786 };
787 
788 struct OpResF_Boundary : public Op_g {
789 
790  OpResF_Boundary(boost::shared_ptr<VectorDouble> &lambda_vals)
791  : Op_g(FVal(), "U", 1), lambdaVals(lambda_vals) {}
792 
793 protected:
794  boost::shared_ptr<VectorDouble> lambdaVals;
795 
796  /**
797  * \brief Integrate local constrains vector
798  */
801  // set size to local vector
802  locVec.resize(nbRows, false);
803  // clear local vector
804  locVec.clear();
805  // get face area
806  const double area = getArea() * bEta;
807  // get integration weight
808  auto t_w = getFTensor0IntegrationWeight();
809  // get base function
810  auto t_u = data.getFTensor0N();
811  // get solution at integration point
812  auto t_lambda = getFTensor0FromVec(*lambdaVals);
813  // make loop over integration points
814  for (int gg = 0; gg != nbIntegrationPts; gg++) {
815  // evaluate function on boundary and scale it by area and integration
816  // weight
817  double alpha = area * t_w;
818  // get element of vector
820  &*locVec.data().begin());
821  for (int rr = 0; rr != nbRows; rr++) {
822  t_a += alpha * t_u * t_lambda;
823  ++t_a;
824  ++t_u;
825  }
826  ++t_w;
827  ++t_lambda;
828  }
830  }
831 };
832 
833 /**
834  * \brief Set integration rule to volume elements
835  *
836  * This rule is used to integrate \f$\nabla v \cdot \nabla u\f$, thus
837  * if approximation field and testing field are polynomial order "p",
838  * then rule for exact integration is 2*(p-1).
839  *
840  * Integration rule is order of polynomial which is calculated exactly. Finite
841  * element selects integration method based on return of this function.
842  *
843  */
844 struct VolRule {
845  int operator()(int, int, int p) const { return 2 * (p - 1); }
846 };
847 
848 /**
849  * \brief Set integration rule to boundary elements
850  *
851  * This is uses to integrate values on the face. Is used to integrate
852  * \f$(\mathbf{n} \cdot \lambda) u\f$, where Lagrange multiplayer
853  * is order "p_row" and approximate function is order "p_col".
854  *
855  * Integration rule is order of polynomial which is calculated exactly. Finite
856  * element selects integration method based on return of this function.
857  *
858  */
859 struct FaceRule {
860  int operator()(int p_row, int p_col, int p_data) const {
861  return 2 * p_data + 1;
862  }
863 };
864 
865 /**
866  * \brief Create finite elements instances
867  *
868  * Create finite element instances and add operators to finite elements.
869  *
870  */
872 
874 
875  /**
876  * \brief Create finite element to calculate matrix and vectors
877  */
879  boost::function<double(const double, const double, const double)> f_u,
880  boost::function<double(const double, const double, const double)>
881  f_source,
882  boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
883  boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
884  boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
885  boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
886  bool trans = true) const {
888 
889  // Create elements element instances
890  domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
892  boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
894  domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
896  boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
898 
899  // Set integration rule to elements instances
900  domain_lhs_fe->getRuleHook = VolRule();
901  domain_rhs_fe->getRuleHook = VolRule();
902  boundary_lhs_fe->getRuleHook = FaceRule();
903  boundary_rhs_fe->getRuleHook = FaceRule();
904 
905  // Add operators to element instances
906  // Add operator grad-grad for calculate matrix
907  domain_lhs_fe->getOpPtrVector().push_back(new OpK());
908  // Add operator to calculate source terms
909  domain_rhs_fe->getOpPtrVector().push_back(new OpF(f_source));
910  // Add operator calculating constrains matrix
911  boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
912  // Add operator calculating constrains vector
913  boundary_rhs_fe->getOpPtrVector().push_back(new Op_g(f_u));
914 
916  }
917 
918  /**
919  * \brief Create finite element to calculate error
920  */
922  boost::function<double(const double, const double, const double)> f_u,
923  boost::function<FTensor::Tensor1<double, 3>(const double, const double,
924  const double)>
925  g_u,
926  Vec global_error,
927  boost::shared_ptr<ForcesAndSourcesCore> &domain_error) const {
929  // Create finite element instance to calculate error
930  domain_error = boost::shared_ptr<ForcesAndSourcesCore>(
932  domain_error->getRuleHook = VolRule();
933  // Set integration rule
934  // Crate shared vector storing values of field "u" on integration points on
935  // element. element is local and is used to exchange data between operators.
936  boost::shared_ptr<VectorDouble> values_at_integration_ptr =
937  boost::make_shared<VectorDouble>();
938  // Storing gradients of field
939  boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
940  boost::make_shared<MatrixDouble>();
941  // Add default operator to calculate field values at integration points
942  domain_error->getOpPtrVector().push_back(
943  new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
944  // Add default operator to calculate field gradient at integration points
945  domain_error->getOpPtrVector().push_back(
946  new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
947  // Add operator to integrate error element by element.
948  domain_error->getOpPtrVector().push_back(
949  new OpError(f_u, g_u, values_at_integration_ptr,
950  grad_at_integration_ptr, global_error));
952  }
953 
954  /**
955  * \brief Create finite element to post-process results
956  */
958  boost::shared_ptr<ForcesAndSourcesCore> &post_proc_volume) const {
959 
961 
962  // Note that user can stack together arbitrary number of operators to
963  // compose complex PDEs.
964 
965  // Post-process results. This is standard element, with functionality
966  // enabling refining mesh for post-processing. In addition in
967  // PostProcOnRefMesh.hpp are implanted set of users operators to
968  // post-processing fields. Here using simplified mechanism for
969  // post-processing finite element, we add operators to save data from field
970  // on mesh tags for ParaView visualization.
971  post_proc_volume = boost::shared_ptr<ForcesAndSourcesCore>(
973  // Add operators to the elements, starting with some generic
974  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
975  post_proc_volume)
976  ->generateReferenceElementMesh();
977  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
978  post_proc_volume)
979  ->addFieldValuesPostProc("U");
980  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
981  post_proc_volume)
982  ->addFieldValuesPostProc("ERROR");
983  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
984  post_proc_volume)
985  ->addFieldValuesGradientPostProc("U");
986 
988  }
989 
990  /**
991  * \brief Create finite element to calculate matrix and vectors
992  */
994  boost::function<double(const double, const double, const double)> f_u,
995  boost::function<double(const double, const double, const double)>
996  f_source,
997  boost::function<double(const double)> a,
998  boost::function<double(const double)> diff_a,
999  boost::shared_ptr<ForcesAndSourcesCore> &domain_lhs_fe,
1000  boost::shared_ptr<ForcesAndSourcesCore> &boundary_lhs_fe,
1001  boost::shared_ptr<ForcesAndSourcesCore> &domain_rhs_fe,
1002  boost::shared_ptr<ForcesAndSourcesCore> &boundary_rhs_fe,
1003  ForcesAndSourcesCore::RuleHookFun vol_rule,
1004  ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule(),
1005  bool trans = true) const {
1007 
1008  // Create elements element instances
1009  domain_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1011  boundary_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1013  domain_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1015  boundary_rhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
1017 
1018  // Set integration rule to elements instances
1019  domain_lhs_fe->getRuleHook = vol_rule;
1020  domain_rhs_fe->getRuleHook = vol_rule;
1021  boundary_lhs_fe->getRuleHook = face_rule;
1022  boundary_rhs_fe->getRuleHook = face_rule;
1023 
1024  // Set integration rule
1025  // Crate shared vector storing values of field "u" on integration points on
1026  // element. element is local and is used to exchange data between operators.
1027  boost::shared_ptr<VectorDouble> values_at_integration_ptr =
1028  boost::make_shared<VectorDouble>();
1029  // Storing gradients of field
1030  boost::shared_ptr<MatrixDouble> grad_at_integration_ptr =
1031  boost::make_shared<MatrixDouble>();
1032  // multipliers values
1033  boost::shared_ptr<VectorDouble> multiplier_at_integration_ptr =
1034  boost::make_shared<VectorDouble>();
1035 
1036  // Add operators to element instances
1037  // Add default operator to calculate field values at integration points
1038  domain_lhs_fe->getOpPtrVector().push_back(
1039  new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
1040  // Add default operator to calculate field gradient at integration points
1041  domain_lhs_fe->getOpPtrVector().push_back(
1042  new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
1043  // Add operator grad-(1+u^2)grad for calculate matrix
1044  domain_lhs_fe->getOpPtrVector().push_back(new OpKt(
1045  a, diff_a, values_at_integration_ptr, grad_at_integration_ptr));
1046 
1047  // Add default operator to calculate field values at integration points
1048  domain_rhs_fe->getOpPtrVector().push_back(
1049  new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
1050  // Add default operator to calculate field gradient at integration points
1051  domain_rhs_fe->getOpPtrVector().push_back(
1052  new OpCalculateScalarFieldGradient<3>("U", grad_at_integration_ptr));
1053  // Add operator to calculate source terms
1054  domain_rhs_fe->getOpPtrVector().push_back(new OpResF_Domain(
1055  f_source, a, values_at_integration_ptr, grad_at_integration_ptr));
1056 
1057  // Add operator calculating constrains matrix
1058  boundary_lhs_fe->getOpPtrVector().push_back(new OpC(trans));
1059 
1060  // Add default operator to calculate field values at integration points
1061  boundary_rhs_fe->getOpPtrVector().push_back(
1062  new OpCalculateScalarFieldValues("U", values_at_integration_ptr));
1063  // Add default operator to calculate values of Lagrange multipliers
1064  boundary_rhs_fe->getOpPtrVector().push_back(
1065  new OpCalculateScalarFieldValues("L", multiplier_at_integration_ptr));
1066  // Add operator calculating constrains vector
1067  boundary_rhs_fe->getOpPtrVector().push_back(
1068  new OpRes_g(f_u, values_at_integration_ptr));
1069  boundary_rhs_fe->getOpPtrVector().push_back(
1070  new OpResF_Boundary(multiplier_at_integration_ptr));
1071 
1073  }
1074 
1075 private:
1077 };
1078 
1079 } // namespace PoissonExample
1080 
1081 #endif //__POISSONOPERATORS_HPP__
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpK
constexpr double alpha
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DataForcesAndSourcesCore::EntData EntData
Deprecated interface functions.
Create finite elements instances.
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< ForcesAndSourcesCore > &post_proc_volume) const
Create finite element to post-process results.
CreateFiniteElements(MoFEM::Interface &m_field)
MoFEMErrorCode createFEToAssembleMatrixAndVectorForNonlinearProblem(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::function< double(const double)> a, boost::function< double(const double)> diff_a, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, ForcesAndSourcesCore::RuleHookFun vol_rule, ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule(), bool trans=true) const
Create finite element to calculate matrix and vectors.
MoFEMErrorCode createFEToEvaluateError(boost::function< double(const double, const double, const double)> f_u, boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> g_u, Vec global_error, boost::shared_ptr< ForcesAndSourcesCore > &domain_error) const
Create finite element to calculate error.
MoFEMErrorCode createFEToAssembleMatrixAndVector(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, bool trans=true) const
Create finite element to calculate matrix and vectors.
Set integration rule to boundary elements.
int operator()(int p_row, int p_col, int p_data) const
Assemble constrains vector.
FTensor::Number< 2 > NZ
z-direction index
Op_g(FVal f_value, const string field_name="L", const double beta=1)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble constrains vectors
boost::function< double(const double, const double, const double)> FVal
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
FTensor::Number< 0 > NX
x-direction index
FTensor::Number< 1 > NY
y-direction index
FVal fValue
Function pointer evaluating values of "U" at the boundary.
template class for integration oh the right hand side
int nbIntegrationPts
number of integration points
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)=0
Class dedicated to assemble operator to global system vector.
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
This function is called by finite element.
OpBaseRhs(const std::string field_name)
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)=0
Class dedicated to integrate operator.
Calculate constrains matrix.
MatrixDouble locMat
local constrains matrix
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
integrate local constrains matrix
OpC(const bool assemble_transpose)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate local constrains matrix.
int nbIntegrationPts
number of integration points
int nbRows
< error code
const bool assembleTranspose
assemble transpose, i.e. CT if set to true
int nbCols
number of dofs on column
Vec globalError
ghost vector with global (integrated over volume) error
UVal uValue
function with exact solution
OpError(UVal u_value, GVal g_value, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals, Vec global_error)
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
boost::function< double(const double, const double, const double)> UVal
FTensor::Number< 2 > NZ
FTensor::Index< 'i', 3 > i
GVal gValue
function with exact solution for gradient
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate error.
FTensor::Number< 1 > NY
boost::function< FTensor::Tensor1< double, 3 > const double, const double, const double)> GVal
FTensor::Number< 0 > NX
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
Assemble error.
Operator calculate source term,.
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
assemble local entity vector to the global right hand side
FTensor::Number< 2 > NZ
OpF(FSource f_source)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double, const double, const double)> FSource
FTensor::Number< 1 > NY
FTensor::Number< 0 > NX
Calculate the grad-grad operator and assemble matrix.
virtual MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
int nbRows
< error code
int nbCols
number if dof on column
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations for give operator.
bool isDiag
true if this block is on diagonal
MatrixDouble locMat
local entity block matrix
int nbIntegrationPts
number of integration points
FTensor::Index< 'i', 3 > i
summit Index
virtual MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble local entity block matrix.
boost::function< double(const double)> diffA
boost::shared_ptr< MatrixDouble > gradVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< VectorDouble > fieldVals
boost::function< double(const double)> A
OpKt(boost::function< double(const double)> a, boost::function< double(const double)> diff_a, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
boost::shared_ptr< VectorDouble > fieldVals
OpRes_g(FVal f_value, boost::shared_ptr< VectorDouble > &field_vals)
boost::shared_ptr< VectorDouble > lambdaVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local constrains vector.
OpResF_Boundary(boost::shared_ptr< VectorDouble > &lambda_vals)
boost::shared_ptr< MatrixDouble > gradVals
boost::shared_ptr< VectorDouble > fieldVals
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Integrate local entity vector.
boost::function< double(const double)> A
FTensor::Index< 'i', 3 > i
OpResF_Domain(FSource f_source, boost::function< double(const double)> a, boost::shared_ptr< VectorDouble > &field_vals, boost::shared_ptr< MatrixDouble > &grad_vals)
Set integration rule to volume elements.
int operator()(int, int, int p) const
Post processing.