v0.13.0
UnsaturatedFlow.hpp
Go to the documentation of this file.
1 /** \file UnsaturatedFlow.hpp
2  * \brief Mix implementation of transport element
3  * \example UnsaturatedFlow.hpp
4  *
5  * \ingroup mofem_mix_transport_elem
6  *
7  */
8 
9 /*
10  * This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #ifndef __UNSATURATD_FLOW_HPP__
25 #define __UNSATURATD_FLOW_HPP__
26 
27 namespace MixTransport {
28 
29 /**
30  * \brief Generic material model for unsaturated water transport
31  *
32  * \note This is abstact class, no physical material model is implemented
33  * here.
34  */
36 
37  virtual ~GenericMaterial() {}
38 
39  static double ePsilon0; ///< Regularization parameter
40  static double ePsilon1; ///< Regularization parameter
41  static double scaleZ; ///< Scale z direction
42 
43  double sCale; ///< Scale time dependent eq.
44 
45  double h; ///< hydraulic head
46  double h_t; ///< rate of hydraulic head
47  // double h_flux_residual; ///< residual at point
48  double K; ///< Hydraulic conductivity [L/s]
49  double diffK; ///< Derivative of hydraulic conductivity [L/s * L^2/F]
50  double C; ///< Capacity [S^2/L^2]
51  double diffC; ///< Derivative of capacity [S^2/L^2 * L^2/F ]
52  double tHeta; ///< Water content
53  double Se; ///< Effective saturation
54 
55  Range tEts; ///< Elements with this material
56 
57  double x, y, z; ///< in meters (L)
58 
59  /**
60  * \brief Initialize head
61  * @return value of head
62  */
63  virtual double initialPcEval() const = 0;
64  virtual void printMatParameters(const int id,
65  const std::string &prefix) const = 0;
66 
67  virtual MoFEMErrorCode calK() {
69  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
70  "Not implemented how to calculate hydraulic conductivity");
72  }
73 
76  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
77  "Not implemented how to calculate derivative of hydraulic "
78  "conductivity");
80  }
81 
82  virtual MoFEMErrorCode calC() {
84  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
85  "Not implemented how to calculate capacity");
87  }
88 
91  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
92  "Not implemented how to calculate capacity");
94  }
95 
98  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
99  "Not implemented how to calculate capacity");
101  }
102 
103  virtual MoFEMErrorCode calSe() {
105  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
106  "Not implemented how to calculate capacity");
108  }
109 };
110 
111 /**
112  * \brief Implementation of operators, problem and finite elements for
113  * unsaturated flow
114  */
116 
117  DM dM; ///< Discrete manager for unsaturated flow problem
118 
120  : MixTransportElement(m_field), dM(PETSC_NULL), lastEvalBcValEnt(0),
122  lastEvalBcBlockFluxId(-1) {}
123 
125  if (dM != PETSC_NULL) {
126  CHKERR DMDestroy(&dM);
127  CHKERRABORT(PETSC_COMM_WORLD, ierr);
128  }
129  }
130 
131  typedef std::map<int, boost::shared_ptr<GenericMaterial>> MaterialsDoubleMap;
132  MaterialsDoubleMap dMatMap; ///< materials database
133 
134  /**
135  * \brief For given element handle get material block Id
136  * @param ent finite element entity handle
137  * @param block_id reference to returned block id
138  * @return error code
139  */
141  int &block_id) const {
143  for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
144  mit != dMatMap.end(); mit++) {
145  if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
146  block_id = mit->first;
148  }
149  }
151  "Element not found, no material data");
153  }
154 
155  /**
156  * \brief Class storing information about boundary condition
157  */
158  struct BcData {
159  Range eNts;
160  double fixValue;
161  boost::function<double(const double x, const double y, const double z)>
163  BcData() : hookFun(NULL) {}
164  };
165  typedef map<int, boost::shared_ptr<BcData>> BcMap;
166  BcMap bcValueMap; ///< Store boundary condition for head capillary pressure
167 
170 
171  /**
172  * \brief Get value on boundary
173  * @param ent entity handle
174  * @param gg number of integration point
175  * @param x x-coordinate
176  * @param y y-coordinate
177  * @param z z-coordinate
178  * @param value returned value
179  * @return error code
180  */
181  MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
182  const double x, const double y, const double z,
183  double &value) {
185  int block_id = -1;
186  if (lastEvalBcValEnt == ent) {
187  block_id = lastEvalBcBlockValId;
188  } else {
189  for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
190  it++) {
191  if (it->second->eNts.find(ent) != it->second->eNts.end()) {
192  block_id = it->first;
193  }
194  }
195  lastEvalBcValEnt = ent;
196  lastEvalBcBlockValId = block_id;
197  }
198  if (block_id >= 0) {
199  if (bcValueMap.at(block_id)->hookFun) {
200  value = bcValueMap.at(block_id)->hookFun(x, y, z);
201  } else {
202  value = bcValueMap.at(block_id)->fixValue;
203  }
204  } else {
205  value = 0;
206  }
208  }
209 
213 
214  /**
215  * \brief essential (Neumann) boundary condition (set fluxes)
216  * @param ent handle to finite element entity
217  * @param x coord
218  * @param y coord
219  * @param z coord
220  * @param flux reference to flux which is set by function
221  * @return [description]
222  */
223  MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
224  const double y, const double z, double &flux) {
226  int block_id = -1;
227  if (lastEvalBcFluxEnt == ent) {
228  block_id = lastEvalBcBlockFluxId;
229  } else {
230  for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
231  it++) {
232  if (it->second->eNts.find(ent) != it->second->eNts.end()) {
233  block_id = it->first;
234  }
235  }
236  lastEvalBcFluxEnt = ent;
237  lastEvalBcBlockFluxId = block_id;
238  }
239  if (block_id >= 0) {
240  if (bcFluxMap.at(block_id)->hookFun) {
241  flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
242  } else {
243  flux = bcFluxMap.at(block_id)->fixValue;
244  }
245  } else {
246  flux = 0;
247  }
249  }
250 
251  /**
252  * \brief Evaluate boundary condition at the boundary
253  */
256 
258  boost::shared_ptr<MethodForForceScaling> valueScale;
259 
260  /**
261  * \brief Constructor
262  */
263  OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name,
264  boost::shared_ptr<MethodForForceScaling> &value_scale)
266  fluxes_name, UserDataOperator::OPROW),
267  cTx(ctx), valueScale(value_scale) {}
268 
269  VectorDouble nF; ///< Vector of residuals
270 
271  /**
272  * \brief Integrate boundary condition
273  * @param side local index of entity
274  * @param type type of entity
275  * @param data data on entity
276  * @return error code
277  */
278  MoFEMErrorCode doWork(int side, EntityType type,
281  if (data.getFieldData().size() == 0)
283  // Get EntityHandle of the finite element
284  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
285  // Resize and clear vector
286  nF.resize(data.getIndices().size());
287  nF.clear();
288  // Get number of integration points
289  int nb_gauss_pts = data.getN().size1();
290  for (int gg = 0; gg < nb_gauss_pts; gg++) {
291  double x, y, z;
292  x = getCoordsAtGaussPts()(gg, 0);
293  y = getCoordsAtGaussPts()(gg, 1);
294  z = getCoordsAtGaussPts()(gg, 2);
295  double value;
296  // get value of boundary condition
297  CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
298  const double w = getGaussPts()(2, gg) * 0.5;
299  const double beta = w * (value - z);
300  noalias(nF) += beta * prod(data.getVectorN<3>(gg), getNormal());
301  }
302  // Scale vector if history evaluating method is given
303  Vec f = getFEMethod()->ts_F;
304  if (valueScale) {
305  CHKERR valueScale->scaleNf(getFEMethod(), nF);
306  }
307  // Assemble vector
308  CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
309  &nF[0], ADD_VALUES);
311  }
312  };
313 
314  /**
315  * \brief Assemble flux residual
316  */
319 
321 
322  OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
324  flux_name, UserDataOperator::OPROW),
325  cTx(ctx) {}
326 
329 
330  MoFEMErrorCode doWork(int side, EntityType type,
333  const int nb_dofs = data.getIndices().size();
334  if (nb_dofs == 0)
336  nF.resize(nb_dofs, false);
337  nF.clear();
338  // Get EntityHandle of the finite element
339  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
340  // Get material block id
341  int block_id;
342  CHKERR cTx.getMaterial(fe_ent, block_id);
343  // Get material block
344  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
345  // block_data->printMatParameters(block_id,"Read material");
346  // Get base function
347  auto t_n_hdiv = data.getFTensor1N<3>();
348  // Get pressure
350  // Get flux
351  auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
352  // Coords at integration points
353  auto t_coords = getFTensor1CoordsAtGaussPts();
354  // Get integration weight
355  auto t_w = getFTensor0IntegrationWeight();
356  // Get volume
357  double vol = getVolume();
358 
360  auto t_base_diff_hdiv = data.getFTensor2DiffN<3, 3>();
361  divVec.resize(nb_dofs, false);
362 
363  // Get material parameters
364  int nb_gauss_pts = data.getN().size1();
365  for (int gg = 0; gg != nb_gauss_pts; gg++) {
366  // Get divergence
367  for(auto &v : divVec) {
368  v = t_base_diff_hdiv(i, i);
369  ++t_base_diff_hdiv;
370  }
371 
372  const double alpha = t_w * vol;
373  block_data->h = t_h;
374  block_data->x = t_coords(0);
375  block_data->y = t_coords(1);
376  block_data->z = t_coords(2);
377  CHKERR block_data->calK();
378  const double K = block_data->K;
379  const double scaleZ = block_data->scaleZ;
380  const double z = t_coords(2); /// z-coordinate at Gauss pt
381  // Calculate pressure gradient
382  noalias(nF) -= alpha * (t_h - z * scaleZ) * divVec;
383  // Calculate presure gradient from flux
384  FTensor::Tensor0<double *> t_nf(&*nF.begin());
385  for (int rr = 0; rr != nb_dofs; rr++) {
386  t_nf += alpha * (1 / K) * (t_n_hdiv(i) * t_flux(i));
387  ++t_n_hdiv; // move to next base function
388  ++t_nf; // move to next element in vector
389  }
390  ++t_h; // move to next integration point
391  ++t_flux;
392  ++t_coords;
393  ++t_w;
394  }
395  // Assemble residual
396  CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs,
397  &*data.getIndices().begin(), &*nF.begin(),
398  ADD_VALUES);
400  }
401  };
402 
403  /**
404  * Assemble mass residual
405  */
408 
410 
411  OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
413  val_name, UserDataOperator::OPROW),
414  cTx(ctx) {}
415 
417 
418  MoFEMErrorCode doWork(int side, EntityType type,
421  const int nb_dofs = data.getIndices().size();
422  if (nb_dofs == 0)
424  // Resize local element vector
425  nF.resize(nb_dofs, false);
426  nF.clear();
427  // Get EntityHandle of the finite element
428  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
429  // Get material block id
430  int block_id;
431  CHKERR cTx.getMaterial(fe_ent, block_id);
432  // Get material block
433  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
434  // Get pressure
436  // Get pressure rate
438  // Flux divergence
439  auto t_div_flux = getFTensor0FromVec(cTx.divergenceAtGaussPts);
440  // Get integration weight
441  auto t_w = getFTensor0IntegrationWeight();
442  // Coords at integration points
443  auto t_coords = getFTensor1CoordsAtGaussPts();
444  // Scale eq.
445  const double scale = block_data->sCale;
446  // Get volume
447  const double vol = getVolume();
448  // Get number of integration points
449  int nb_gauss_pts = data.getN().size1();
450  for (int gg = 0; gg != nb_gauss_pts; gg++) {
451  const double alpha = t_w * vol * scale;
452  block_data->h = t_h;
453  block_data->x = t_coords(0);
454  block_data->y = t_coords(1);
455  block_data->z = t_coords(2);
456  CHKERR block_data->calC();
457  const double C = block_data->C;
458  // Calculate flux conservation
459  noalias(nF) += (alpha * (t_div_flux + C * t_h_t)) * data.getN(gg);
460  ++t_h;
461  ++t_h_t;
462  ++t_div_flux;
463  ++t_coords;
464  ++t_w;
465  }
466  // Assemble local vector
467  Vec f = getFEMethod()->ts_F;
468  CHKERR VecSetValues(f, nb_dofs, &*data.getIndices().begin(), &*nF.begin(),
469  ADD_VALUES);
471  }
472  };
473 
476 
478 
480  const std::string flux_name)
482  flux_name, flux_name, UserDataOperator::OPROWCOL),
483  cTx(ctx) {
484  sYmm = true;
485  }
486 
488 
490 
491  /**
492  * \brief Assemble matrix
493  * @param row_side local index of row entity on element
494  * @param col_side local index of col entity on element
495  * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
496  * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
497  * @param row_data data for row
498  * @param col_data data for col
499  * @return error code
500  */
501  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
502  EntityType col_type,
503  EntitiesFieldData::EntData &row_data,
504  EntitiesFieldData::EntData &col_data) {
506  const int nb_row = row_data.getIndices().size();
507  const int nb_col = col_data.getIndices().size();
508  if (nb_row == 0)
510  if (nb_col == 0)
512  nN.resize(nb_row, nb_col, false);
513  nN.clear();
514  // Get EntityHandle of the finite element
515  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
516  // Get material block id
517  int block_id;
518  CHKERR cTx.getMaterial(fe_ent, block_id);
519  // Get material block
520  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
521  // Get pressure
523  // Coords at integration points
524  auto t_coords = getFTensor1CoordsAtGaussPts();
525  // Get base functions
526  auto t_n_hdiv_row = row_data.getFTensor1N<3>();
527  // Get integration weight
528  auto t_w = getFTensor0IntegrationWeight();
529  // Get volume
530  const double vol = getVolume();
531  int nb_gauss_pts = row_data.getN().size1();
532  for (int gg = 0; gg != nb_gauss_pts; gg++) {
533  block_data->h = t_h;
534  block_data->x = t_coords(0);
535  block_data->y = t_coords(1);
536  block_data->z = t_coords(2);
537  CHKERR block_data->calK();
538  const double K = block_data->K;
539  // get integration weight and multiply by element volume
540  const double alpha = t_w * vol;
541  const double beta = alpha * (1 / K);
542  FTensor::Tensor0<double *> t_a(&*nN.data().begin());
543  for (int kk = 0; kk != nb_row; kk++) {
544  auto t_n_hdiv_col = col_data.getFTensor1N<3>(gg, 0);
545  for (int ll = 0; ll != nb_col; ll++) {
546  t_a += beta * (t_n_hdiv_row(j) * t_n_hdiv_col(j));
547  ++t_n_hdiv_col;
548  ++t_a;
549  }
550  ++t_n_hdiv_row;
551  }
552  ++t_coords;
553  ++t_h;
554  ++t_w;
555  }
556  Mat a = getFEMethod()->ts_B;
557  CHKERR MatSetValues(a, nb_row, &*row_data.getIndices().begin(), nb_col,
558  &*col_data.getIndices().begin(), &*nN.data().begin(),
559  ADD_VALUES);
560  // matrix is symmetric, assemble other part
561  if (row_side != col_side || row_type != col_type) {
562  transNN.resize(nb_col, nb_row);
563  noalias(transNN) = trans(nN);
564  CHKERR MatSetValues(a, nb_col, &*col_data.getIndices().begin(), nb_row,
565  &*row_data.getIndices().begin(),
566  &*transNN.data().begin(), ADD_VALUES);
567  }
569  }
570  };
571 
572  struct OpVU_L2L2
574 
576 
577  OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
579  value_name, value_name, UserDataOperator::OPROWCOL),
580  cTx(ctx) {
581  sYmm = true;
582  }
583 
585 
586  /**
587  * \brief Assemble matrix
588  * @param row_side local index of row entity on element
589  * @param col_side local index of col entity on element
590  * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
591  * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
592  * @param row_data data for row
593  * @param col_data data for col
594  * @return error code
595  */
596  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
597  EntityType col_type,
598  EntitiesFieldData::EntData &row_data,
599  EntitiesFieldData::EntData &col_data) {
601  int nb_row = row_data.getIndices().size();
602  int nb_col = col_data.getIndices().size();
603  if (nb_row == 0)
605  if (nb_col == 0)
607  nN.resize(nb_row, nb_col, false);
608  nN.clear();
609  // Get EntityHandle of the finite element
610  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
611  // Get material block id
612  int block_id;
613  CHKERR cTx.getMaterial(fe_ent, block_id);
614  // Get material block
615  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
616  // Get pressure
618  // // Get pressure
619  // auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
620  // Get pressure rate
622  // Get integration weight
623  auto t_w = getFTensor0IntegrationWeight();
624  // Coords at integration points
625  auto t_coords = getFTensor1CoordsAtGaussPts();
626  // Scale eq.
627  const double scale = block_data->sCale;
628  // Time step factor
629  double ts_a = getFEMethod()->ts_a;
630  // get volume
631  const double vol = getVolume();
632  int nb_gauss_pts = row_data.getN().size1();
633  // get base functions on rows
634  auto t_n_row = row_data.getFTensor0N();
635  for (int gg = 0; gg != nb_gauss_pts; gg++) {
636  // get integration weight and multiply by element volume
637  double alpha = t_w * vol * scale;
638  // evaluate material model at integration points
639  // to calculate capacity and tangent of capacity term
640  block_data->h = t_h;
641  block_data->h_t = t_h_t;
642  block_data->x = t_coords(0);
643  block_data->y = t_coords(1);
644  block_data->z = t_coords(2);
645  CHKERR block_data->calC();
646  CHKERR block_data->calDiffC();
647  const double C = block_data->C;
648  const double diffC = block_data->diffC;
649  // assemble local entity tangent matrix
651  &*nN.data().begin());
652  // iterate base functions on rows
653  for (int kk = 0; kk != nb_row; kk++) {
654  // get first base function on column at integration point gg
655  auto t_n_col = col_data.getFTensor0N(gg, 0);
656  // iterate base functions on columns
657  for (int ll = 0; ll != nb_col; ll++) {
658  // assemble elements of local matrix
659  t_a += (alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
660  ++t_n_col; // move to next base function on column
661  ++t_a; // move to next element in local tangent matrix
662  }
663  ++t_n_row; // move to next base function on row
664  }
665  ++t_w; // move to next integration weight
666  ++t_coords; // move to next coordinate at integration point
667  ++t_h; // move to next capillary head at integration point
668  // ++t_flux_residual;
669  ++t_h_t; // move to next capillary head rate at integration point
670  }
671  Mat a = getFEMethod()->ts_B;
672  CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
673  &col_data.getIndices()[0], &*nN.data().begin(),
674  ADD_VALUES);
676  }
677  };
678 
681 
683 
684  /**
685  * \brief Constructor
686  */
688  const std::string &val_name_row,
689  const std::string &flux_name_col)
691  val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
692  cTx(ctx) {}
693 
696 
697  /**
698  * \brief Do calculations
699  * @param row_side local index of entity on row
700  * @param col_side local index of entity on column
701  * @param row_type type of row entity
702  * @param col_type type of col entity
703  * @param row_data row data structure carrying information about base
704  * functions, DOFs indices, etc.
705  * @param col_data column data structure carrying information about base
706  * functions, DOFs indices, etc.
707  * @return error code
708  */
709  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
710  EntityType col_type,
711  EntitiesFieldData::EntData &row_data,
712  EntitiesFieldData::EntData &col_data) {
714  int nb_row = row_data.getFieldData().size();
715  int nb_col = col_data.getFieldData().size();
716  if (nb_row == 0)
718  if (nb_col == 0)
720  // Get EntityHandle of the finite element
721  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
722  // Get material block id
723  int block_id;
724  CHKERR cTx.getMaterial(fe_ent, block_id);
725  // Get material block
726  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
727  nN.resize(nb_row, nb_col, false);
728  divVec.resize(nb_col, false);
729  nN.clear();
730  // take derivatives of base functions
732  auto t_base_diff_hdiv = col_data.getFTensor2DiffN<3, 3>();
733  // Scale eq.
734  const double scale = block_data->sCale;
735  int nb_gauss_pts = row_data.getN().size1();
736  for (int gg = 0; gg < nb_gauss_pts; gg++) {
737  double alpha = getGaussPts()(3, gg) * getVolume() * scale;
738  for(auto &v : divVec) {
739  v = t_base_diff_hdiv(i, i);
740  ++t_base_diff_hdiv;
741  }
742  noalias(nN) += alpha * outer_prod(row_data.getN(gg), divVec);
743  }
744  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
745  &row_data.getIndices()[0], nb_col,
746  &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
748  }
749  };
750 
753 
755 
756  /**
757  * \brief Constructor
758  */
760  const std::string &flux_name_col,
761  const std::string &val_name_row)
763  flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
764  cTx(ctx) {}
765 
769 
770  /**
771  * \brief Do calculations
772  * @param row_side local index of entity on row
773  * @param col_side local index of entity on column
774  * @param row_type type of row entity
775  * @param col_type type of col entity
776  * @param row_data row data structure carrying information about base
777  * functions, DOFs indices, etc.
778  * @param col_data column data structure carrying information about base
779  * functions, DOFs indices, etc.
780  * @return error code
781  */
782  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
783  EntityType col_type,
784  EntitiesFieldData::EntData &row_data,
785  EntitiesFieldData::EntData &col_data) {
787  int nb_row = row_data.getFieldData().size();
788  int nb_col = col_data.getFieldData().size();
789  if (nb_row == 0)
791  if (nb_col == 0)
793  nN.resize(nb_row, nb_col, false);
794  divVec.resize(nb_row, false);
795  nN.clear();
796  // Get EntityHandle of the finite element
797  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
798  // Get material block id
799  int block_id;
800  CHKERR cTx.getMaterial(fe_ent, block_id);
801  // Get material block
802  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
803  // Get pressure
805  // Get flux
806  auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
807  // Coords at integration points
808  auto t_coords = getFTensor1CoordsAtGaussPts();
809  // Get integration weight
810  auto t_w = getFTensor0IntegrationWeight();
811  // Get base function
812  auto t_n_hdiv_row = row_data.getFTensor1N<3>();
813  // Get derivatives of base functions
815  auto t_base_diff_hdiv = row_data.getFTensor2DiffN<3, 3>();
816  // Get volume
817  double vol = getVolume();
818  int nb_gauss_pts = row_data.getN().size1();
819  for (int gg = 0; gg != nb_gauss_pts; gg++) {
820  block_data->h = t_h;
821  block_data->x = t_coords(0);
822  block_data->y = t_coords(1);
823  block_data->z = t_coords(2);
824  CHKERR block_data->calK();
825  CHKERR block_data->calDiffK();
826  const double K = block_data->K;
827  // const double z = t_coords(2);
828  const double KK = K * K;
829  const double diffK = block_data->diffK;
830  double alpha = t_w * vol;
831  for(auto &v : divVec) {
832  v = t_base_diff_hdiv(i, i);
833  ++t_base_diff_hdiv;
834  }
835  noalias(nN) -= alpha * outer_prod(divVec, col_data.getN(gg));
836  FTensor::Tensor0<double *> t_a(&*nN.data().begin());
837  for (int rr = 0; rr != nb_row; rr++) {
838  double beta = alpha * (-diffK / KK) * (t_n_hdiv_row(i) * t_flux(i));
839  auto t_n_col = col_data.getFTensor0N(gg, 0);
840  for (int cc = 0; cc != nb_col; cc++) {
841  t_a += beta * t_n_col;
842  ++t_n_col;
843  ++t_a;
844  }
845  ++t_n_hdiv_row;
846  }
847  ++t_w;
848  ++t_coords;
849  ++t_h;
850  ++t_flux;
851  }
852  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
853  &row_data.getIndices()[0], nb_col,
854  &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
856  }
857  };
858 
863  const std::string &val_name)
865  val_name, UserDataOperator::OPROW),
866  cTx(ctx) {}
867 
870 
871  MoFEMErrorCode doWork(int side, EntityType type,
874  if (data.getFieldData().size() == 0)
876  int nb_dofs = data.getFieldData().size();
877  int nb_gauss_pts = data.getN().size1();
878  if (nb_dofs != static_cast<int>(data.getN().size2())) {
879  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
880  "wrong number of dofs");
881  }
882  nN.resize(nb_dofs, nb_dofs, false);
883  nF.resize(nb_dofs, false);
884  nN.clear();
885  nF.clear();
886 
887  // Get EntityHandle of the finite element
888  EntityHandle fe_ent = getFEEntityHandle();
889  // Get material block id
890  int block_id;
891  CHKERR cTx.getMaterial(fe_ent, block_id);
892  // Get material block
893  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
894 
895  // loop over integration points
896  for (int gg = 0; gg < nb_gauss_pts; gg++) {
897  // get coordinates at integration point
898  block_data->x = getCoordsAtGaussPts()(gg, 0);
899  block_data->y = getCoordsAtGaussPts()(gg, 1);
900  block_data->z = getCoordsAtGaussPts()(gg, 2);
901  // get weight for integration rule
902  double alpha = getGaussPts()(2, gg) * getVolume();
903  nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
904  nF += alpha * block_data->initialPcEval() * data.getN(gg);
905  }
906 
907  // factor matrix
909  // solve local problem
910  cholesky_solve(nN, nF, ublas::lower());
911 
912  // set solution to vector
913  CHKERR VecSetValues(cTx.D1, nb_dofs, &*data.getIndices().begin(),
914  &*nF.begin(), INSERT_VALUES);
915 
917  }
918  };
919 
923 
924  /**
925  * \brief Constructor
926  */
928  const std::string fluxes_name)
930  fluxes_name, UserDataOperator::OPROW),
931  cTx(ctx) {}
932 
934 
935  /**
936  * \brief Integrate boundary condition
937  * @param side local index of entity
938  * @param type type of entity
939  * @param data data on entity
940  * @return error code
941  */
942  MoFEMErrorCode doWork(int side, EntityType type,
945  int nb_dofs = data.getFieldData().size();
946  if (nb_dofs == 0)
948  // Get base function
949  auto t_n_hdiv = data.getFTensor1N<3>();
950  // get normal of face
951  auto t_normal = getFTensor1Normal();
952  // Integration weight
953  auto t_w = getFTensor0IntegrationWeight();
954  double flux_on_entity = 0;
955  int nb_gauss_pts = data.getN().size1();
956  for (int gg = 0; gg < nb_gauss_pts; gg++) {
957  auto t_data = data.getFTensor0FieldData();
958  for (int rr = 0; rr != nb_dofs; rr++) {
959  flux_on_entity -= (0.5 * t_data * t_w) * (t_n_hdiv(i) * t_normal(i));
960  ++t_n_hdiv;
961  ++t_data;
962  }
963  ++t_w;
964  }
965  CHKERR VecSetValue(cTx.ghostFlux, 0, flux_on_entity, ADD_VALUES);
967  }
968  };
969 
970  /**
971  * Operator used to post-process results for unsaturated infiltration problem.
972  * Operator should with element for post-processing results, i.e.
973  * PostProcVolumeOnRefinedMesh
974  */
977 
980  std::vector<EntityHandle> &mapGaussPts;
981 
983  moab::Interface &post_proc_mesh,
984  std::vector<EntityHandle> &map_gauss_pts,
985  const std::string field_name)
988  cTx(ctx), postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
989 
990  MoFEMErrorCode doWork(int side, EntityType type,
993  int nb_dofs = data.getFieldData().size();
994  if (nb_dofs == 0)
996 
997  // Get EntityHandle of the finite element
998  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
999  // Get material block id
1000  int block_id;
1001  CHKERR cTx.getMaterial(fe_ent, block_id);
1002  // Get material block
1003  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
1004 
1005  // Set bloc Id
1006  Tag th_id;
1007  int def_block_id = -1;
1008  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
1009  MB_TAG_CREAT | MB_TAG_SPARSE,
1010  &def_block_id);
1011 
1012  // Create mesh tag. Tags are created on post-processing mesh and
1013  // visable in post-processor, e.g. Paraview
1014  double zero = 0;
1015  Tag th_theta;
1016  CHKERR postProcMesh.tag_get_handle("THETA", 1, MB_TYPE_DOUBLE, th_theta,
1017  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1018  Tag th_se;
1019  CHKERR postProcMesh.tag_get_handle("Se", 1, MB_TYPE_DOUBLE, th_se,
1020  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1021  // Tag th_ks;
1022  // CHKERR postProcMesh.tag_get_handle(
1023  // "Ks",1,MB_TYPE_DOUBLE,th_ks,
1024  // MB_TAG_CREAT|MB_TAG_SPARSE,&zero
1025  // );
1026  // CHKERR postProcMesh.tag_set_data(th_ks,&fe_ent,1,&block_data->Ks);
1027  Tag th_k;
1028  CHKERR postProcMesh.tag_get_handle("K", 1, MB_TYPE_DOUBLE, th_k,
1029  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1030  Tag th_c;
1031  CHKERR postProcMesh.tag_get_handle("C", 1, MB_TYPE_DOUBLE, th_c,
1032  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1033 
1034  // Get pressure at integration points
1036  // Coords at integration points
1037  auto t_coords = getFTensor1CoordsAtGaussPts();
1038 
1039  int nb_gauss_pts = data.getN().size1();
1040  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1041  block_data->h = t_h;
1042  block_data->x = t_coords(0);
1043  block_data->y = t_coords(1);
1044  block_data->z = t_coords(2);
1045  // Calculate theta (water content) and save it on mesh tags
1046  CHKERR block_data->calTheta();
1047  double theta = block_data->tHeta;
1048  CHKERR postProcMesh.tag_set_data(th_theta, &mapGaussPts[gg], 1, &theta);
1049  CHKERR block_data->calSe();
1050  // Calculate Se (effective saturation and save it on the mesh tags)
1051  double Se = block_data->Se;
1052  CHKERR postProcMesh.tag_set_data(th_se, &mapGaussPts[gg], 1, &Se);
1053  // Calculate K (hydraulic conductivity) and save it on the mesh tags
1054  CHKERR block_data->calK();
1055  double K = block_data->K;
1056  CHKERR postProcMesh.tag_set_data(th_k, &mapGaussPts[gg], 1, &K);
1057  // Calculate water capacity and save it on the mesh tags
1058  CHKERR block_data->calC();
1059  double C = block_data->C;
1060  CHKERR postProcMesh.tag_set_data(th_c, &mapGaussPts[gg], 1, &C);
1061 
1062  // Set block iD
1063  CHKERR postProcMesh.tag_set_data(th_id, &mapGaussPts[gg], 1, &block_id);
1064 
1065  ++t_h;
1066  ++t_coords;
1067  }
1068 
1070  }
1071  };
1072 
1073  /**
1074  * Finite element implementation called by TS monitor. Element calls other
1075  * finite elements to evaluate material properties and save results on the
1076  * mesh.
1077  *
1078  * \note Element overloaded only FEMethod::postProcess methos where other
1079  * elements are called.
1080  */
1082 
1084  boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProc;
1085  boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
1086 
1087  const int fRequency;
1088 
1090  boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc,
1091  boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
1092  const int frequency)
1093  : cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
1094  fRequency(frequency) {}
1095 
1099  }
1100 
1104  }
1105 
1108 
1109  // Get time step
1110  int step;
1111  CHKERR TSGetTimeStepNumber(ts, &step);
1112 
1113  if ((step) % fRequency == 0) {
1114  // Post-process results and save in the file
1115  PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
1116  fRequency);
1118  CHKERR postProc->writeFile(
1119  string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
1120  }
1121 
1122  // Integrate fluxes on faces where pressure head is applied
1123  CHKERR VecZeroEntries(cTx.ghostFlux);
1124  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1125  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1126  // Run finite element to integrate fluxes
1128  CHKERR VecAssemblyBegin(cTx.ghostFlux);
1129  CHKERR VecAssemblyEnd(cTx.ghostFlux);
1130  // accumulate errors from processors
1131  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1132  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1133  // scatter errors to all processors
1134  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1135  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1136  double *ghost_flux;
1137  CHKERR VecGetArray(cTx.ghostFlux, &ghost_flux);
1138  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Flux at time %6.4g %6.4g\n", ts_t,
1139  ghost_flux[0]);
1140  CHKERR VecRestoreArray(cTx.ghostFlux, &ghost_flux);
1141 
1143  }
1144  };
1145 
1146  /// \brief add fields
1147  MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
1148  const int order) {
1150  // Fields
1153  CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
1154  // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
1155 
1156  // meshset consisting all entities in mesh
1157  EntityHandle root_set = mField.get_moab().get_root_set();
1158  // add entities to field
1159 
1161  if (it->getName().compare(0, 4, "SOIL") != 0)
1162  continue;
1163  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1164  MBTET, fluxes);
1165  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1166  MBTET, values);
1167  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1168  MBTET, values + "_t");
1169  // CHKERR mField.add_ents_to_field_by_type(
1170  // dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
1171  // );
1172  }
1173 
1174  CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
1175  CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
1176  CHKERR mField.set_field_order(root_set, MBTET, values, order);
1177  CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
1178  // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
1180  }
1181 
1182  /// \brief add finite elements
1183  MoFEMErrorCode addFiniteElements(const std::string &fluxes_name,
1184  const std::string &values_name) {
1186 
1187  // Define element "MIX". Note that this element will work with fluxes_name
1188  // and values_name. This reflect bilinear form for the problem
1197  values_name + "_t");
1198  // CHKERR
1199  // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
1200 
1202  if (it->getName().compare(0, 4, "SOIL") != 0)
1203  continue;
1205  dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
1206  }
1207 
1208  // Define element to integrate natural boundary conditions, i.e. set values.
1209  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
1211  fluxes_name);
1213  fluxes_name);
1215  fluxes_name);
1217  values_name);
1218 
1220  if (it->getName().compare(0, 4, "HEAD") != 0)
1221  continue;
1223  bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
1224  }
1225 
1226  // Define element to apply essential boundary conditions.
1227  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
1229  fluxes_name);
1231  fluxes_name);
1233  fluxes_name);
1235  values_name);
1236 
1238  if (it->getName().compare(0, 4, "FLUX") != 0)
1239  continue;
1241  bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
1242  }
1243 
1245  }
1246 
1247  /**
1248  * \brief Build problem
1249  * @param ref_level mesh refinement on which mesh problem you like to built.
1250  * @return error code
1251  */
1252  MoFEMErrorCode buildProblem(Range zero_flux_ents,
1253  BitRefLevel ref_level = BitRefLevel().set(0)) {
1255 
1256  // Build fields
1258  // Build finite elements
1260  CHKERR mField.build_finite_elements("MIX_BCFLUX");
1261  CHKERR mField.build_finite_elements("MIX_BCVALUE");
1262  // Build adjacencies of degrees of freedom and elements
1263  CHKERR mField.build_adjacencies(ref_level);
1264 
1265  // create DM instance
1266  CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
1267  // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
1268  CHKERR DMSetType(dM, "DMMOFEM");
1269  // mesh is portioned, each process keeps only part of problem
1270  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1271  // creates problem in DM
1272  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
1273  // discretised problem creates square matrix (that makes some optimizations)
1274  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1275  // set DM options from command line
1276  CHKERR DMSetFromOptions(dM);
1277  // add finite elements
1278  CHKERR DMMoFEMAddElement(dM, "MIX");
1279  CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
1280  CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
1281  // constructor data structures
1282  CHKERR DMSetUp(dM);
1283 
1284  // remove zero flux dofs
1285  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
1286  "MIX", "FLUXES", zero_flux_ents);
1287 
1288  PetscSection section;
1289  CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
1290  CHKERR DMSetSection(dM, section);
1291  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
1292  CHKERR PetscSectionDestroy(&section);
1293 
1295  }
1296 
1297  boost::shared_ptr<ForcesAndSourcesCore>
1298  feFaceBc; ///< Elemnet to calculate essential bc
1299  boost::shared_ptr<ForcesAndSourcesCore>
1300  feFaceRhs; ///< Face element apply natural bc
1301  boost::shared_ptr<ForcesAndSourcesCore>
1302  feVolInitialPc; ///< Calculate inital boundary conditions
1303  boost::shared_ptr<ForcesAndSourcesCore>
1304  feVolRhs; ///< Assemble residual vector
1305  boost::shared_ptr<ForcesAndSourcesCore> feVolLhs; ///< Assemble tangent matrix
1306  boost::shared_ptr<MethodForForceScaling>
1307  scaleMethodFlux; ///< Method scaling fluxes
1308  boost::shared_ptr<MethodForForceScaling>
1309  scaleMethodValue; ///< Method scaling values
1310  boost::shared_ptr<FEMethod> tsMonitor; ///< Element used by TS monitor to
1311  ///< postprocess results at time step
1312 
1313  boost::shared_ptr<VectorDouble>
1314  headRateAtGaussPts; ///< Vector keeps head rate
1315  // boost::shared_ptr<VectorDouble> resAtGaussPts; ///< Residual field
1316 
1317  /**
1318  * \brief Set integration rule to volume elements
1319  *
1320  */
1321  struct VolRule {
1322  int operator()(int, int, int p_data) const { return 2 * p_data + p_data; }
1323  };
1324  /**
1325  * \brief Set integration rule to boundary elements
1326  *
1327  */
1328  struct FaceRule {
1329  int operator()(int p_row, int p_col, int p_data) const {
1330  return 2 * p_data;
1331  }
1332  };
1333 
1334  std::vector<int> bcVecIds;
1336 
1337  /**
1338  * \brief Pre-peprocessing
1339  * Set head pressute rate and get inital essential boundary conditions
1340  */
1341  struct preProcessVol {
1343  boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1344  // std::string mArk;
1345 
1348  boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
1349  )
1350  : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1353  // Update pressure rates
1354  CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1355  fePtr->problemPtr, "VALUES", string("VALUES") + "_t", ROW,
1356  fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1357  switch (fePtr->ts_ctx) {
1358  case TSMethod::CTX_TSSETIFUNCTION:
1359  CHKERR VecSetOption(fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
1360  PETSC_TRUE);
1361  if (!cTx.bcIndices.empty()) {
1362  double scale;
1363  CHKERR cTx.scaleMethodFlux->getForceScale(fePtr->ts_t, scale);
1364  if (cTx.bcVecIds.size() != cTx.bcIndices.size()) {
1365  cTx.bcVecIds.insert(cTx.bcVecIds.begin(), cTx.bcIndices.begin(),
1366  cTx.bcIndices.end());
1367  cTx.bcVecVals.resize(cTx.bcVecIds.size(), false);
1368  cTx.vecValsOnBc.resize(cTx.bcVecIds.size(), false);
1369  }
1370  CHKERR VecGetValues(cTx.D0, cTx.bcVecIds.size(),
1371  &*cTx.bcVecIds.begin(), &*cTx.bcVecVals.begin());
1372  CHKERR VecGetValues(fePtr->ts_u, cTx.bcVecIds.size(),
1373  &*cTx.bcVecIds.begin(),
1374  &*cTx.vecValsOnBc.begin());
1375  cTx.bcVecVals *= scale;
1376  VectorDouble::iterator vit = cTx.bcVecVals.begin();
1377  const NumeredDofEntity *dof_ptr;
1378  for (std::vector<int>::iterator it = cTx.bcVecIds.begin();
1379  it != cTx.bcVecIds.end(); it++, vit++) {
1380  if (auto dof_ptr =
1381  fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
1382  .lock())
1383  dof_ptr->getFieldData() = *vit;
1384  }
1385  } else {
1386  cTx.bcVecIds.resize(0);
1387  cTx.bcVecVals.resize(0);
1388  cTx.vecValsOnBc.resize(0);
1389  }
1390  break;
1391  default:
1392  // don nothing
1393  break;
1394  }
1396  }
1397  };
1398 
1399  /**
1400  * \brief Post proces method for volume element
1401  * Assemble vectors and matrices and apply essential boundary conditions
1402  */
1405  boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1406  // std::string mArk;
1409  boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
1410  )
1411  : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1414  switch (fePtr->ts_ctx) {
1415  case TSMethod::CTX_TSSETIJACOBIAN: {
1416  CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1417  CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1418  // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1419  // std::string wait;
1420  // std::cin >> wait;
1421  CHKERR MatZeroRowsColumns(fePtr->ts_B, cTx.bcVecIds.size(),
1422  &*cTx.bcVecIds.begin(), 1, PETSC_NULL,
1423  PETSC_NULL);
1424  CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1425  CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1426  // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1427  // std::string wait;
1428  // std::cin >> wait;
1429  } break;
1430  case TSMethod::CTX_TSSETIFUNCTION: {
1431  CHKERR VecAssemblyBegin(fePtr->ts_F);
1432  CHKERR VecAssemblyEnd(fePtr->ts_F);
1433  if (!cTx.bcVecIds.empty()) {
1435  CHKERR VecSetValues(fePtr->ts_F, cTx.bcVecIds.size(),
1436  &*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
1437  INSERT_VALUES);
1438  }
1439  CHKERR VecAssemblyBegin(fePtr->ts_F);
1440  CHKERR VecAssemblyEnd(fePtr->ts_F);
1441  } break;
1442  default:
1443  // don nothing
1444  break;
1445  }
1447  }
1448  };
1449 
1450  /**
1451  * \brief Create finite element instances
1452  * @param vol_rule integration rule for volume element
1453  * @param face_rule integration rule for boundary element
1454  * @return error code
1455  */
1457  setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule = VolRule(),
1458  ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule()) {
1460 
1461  // create finite element instances
1462  feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
1464  feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1466  feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
1468  feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
1470  feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1472  // set integration rule to elements
1473  feFaceBc->getRuleHook = face_rule;
1474  feFaceRhs->getRuleHook = face_rule;
1475  feVolInitialPc->getRuleHook = vol_rule;
1476  feVolLhs->getRuleHook = vol_rule;
1477  feVolRhs->getRuleHook = vol_rule;
1478  // set function hook for finite element postprocessing stage
1479  feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
1480  feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
1481  feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
1482  feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
1483 
1484  // Set Piola transform
1485  feFaceBc->getOpPtrVector().push_back(
1486  new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1487  feFaceRhs->getOpPtrVector().push_back(
1488  new OpHOSetContravariantPiolaTransformOnFace3D(HDIV));
1489 
1490  // create method for setting history for fluxes on boundary
1491  scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
1492  new TimeForceScale("-flux_history", false));
1493 
1494  // create method for setting history for presure heads on boundary
1495  scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
1496  new TimeForceScale("-head_history", false));
1497 
1498 
1499  // Set operator to calculate essential boundary conditions
1500  feFaceBc->getOpPtrVector().push_back(
1501  new OpEvaluateBcOnFluxes(*this, "FLUXES"));
1502 
1503  // Set operator to calculate initial capillary pressure
1504  feVolInitialPc->getOpPtrVector().push_back(
1505  new OpEvaluateInitiallHead(*this, "VALUES"));
1506 
1507  // set residual face from Neumann terms, i.e. applied pressure
1508  feFaceRhs->getOpPtrVector().push_back(
1509  new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
1510  // set residual finite element operators
1511  headRateAtGaussPts = boost::make_shared<VectorDouble>();
1512  // resAtGaussPts = boost::make_shared<VectorDouble>();
1513  feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1514  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1515  feVolRhs->getOpPtrVector().push_back(
1516  new OpValuesAtGaussPts(*this, "VALUES"));
1517  feVolRhs->getOpPtrVector().push_back(
1518  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1519  feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
1520  feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
1521  feVolRhs->getOpPtrVector().back().opType =
1522  ForcesAndSourcesCore::UserDataOperator::OPROW;
1523  // set tangent matrix finite element operators
1524  feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1525  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1526  // feVolLhs->getOpPtrVector().push_back(
1527  // new
1528  // OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
1529  // );
1530  feVolLhs->getOpPtrVector().push_back(
1531  new OpValuesAtGaussPts(*this, "VALUES"));
1532  feVolLhs->getOpPtrVector().push_back(
1533  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1534  feVolLhs->getOpPtrVector().push_back(
1535  new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
1536  feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
1537  feVolLhs->getOpPtrVector().push_back(
1538  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
1539  feVolLhs->getOpPtrVector().push_back(
1540  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
1541 
1542  // Adding finite elements to DM, time solver will ask for it to assemble
1543  // tangent matrices and residuals
1544  boost::shared_ptr<FEMethod> null;
1545  CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
1546  CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
1547  CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
1548 
1549  // setting up post-processing
1550  boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_process(
1552  CHKERR post_process->generateReferenceElementMesh();
1553  CHKERR post_process->addFieldValuesPostProc("VALUES");
1554  CHKERR post_process->addFieldValuesPostProc("VALUES_t");
1555  // CHKERR post_process->addFieldValuesPostProc("FLUXES_residual");
1556  CHKERR post_process->addFieldValuesPostProc("FLUXES");
1557  post_process->getOpPtrVector().push_back(
1558  new OpValuesAtGaussPts(*this, "VALUES"));
1559  post_process->getOpPtrVector().push_back(
1560  new OpPostProcMaterial(*this, post_process->postProcMesh,
1561  post_process->mapGaussPts, "VALUES"));
1562 
1563  // Integrate fluxes on boundary
1564  boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
1565  flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
1567  flux_integrate->getOpPtrVector().push_back(
1568  new OpIntegrateFluxes(*this, "FLUXES"));
1569  int frequency = 1;
1570  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
1571  CHKERR PetscOptionsInt(
1572  "-how_often_output",
1573  "frequency how often results are dumped on hard disk", "", frequency,
1574  &frequency, NULL);
1575  ierr = PetscOptionsEnd();
1576  CHKERRG(ierr);
1577 
1578  tsMonitor = boost::shared_ptr<FEMethod>(
1579  new MonitorPostProc(*this, post_process, flux_integrate, frequency));
1580  TsCtx *ts_ctx;
1581  DMMoFEMGetTsCtx(dM, &ts_ctx);
1582  ts_ctx->getPostProcessMonitor().push_back(tsMonitor);
1584  }
1585 
1586  Vec D1; ///< Vector with inital head capillary pressure
1587  Vec ghostFlux; ///< Ghost Vector of integrated fluxes
1588 
1589  /**
1590  * \brief Create vectors and matrices
1591  * @return Error code
1592  */
1595  CHKERR DMCreateMatrix(dM, &Aij);
1596  CHKERR DMCreateGlobalVector(dM, &D0);
1597  CHKERR VecDuplicate(D0, &D1);
1598  CHKERR VecDuplicate(D0, &D);
1599  CHKERR VecDuplicate(D0, &F);
1600  int ghosts[] = {0};
1601  int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
1602  int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
1603  CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
1604  &ghostFlux);
1606  }
1607 
1608  /**
1609  * \brief Delete matrices and vector when no longer needed
1610  * @return error code
1611  */
1614  CHKERR MatDestroy(&Aij);
1615  CHKERR VecDestroy(&D);
1616  CHKERR VecDestroy(&D0);
1617  CHKERR VecDestroy(&D1);
1618  CHKERR VecDestroy(&F);
1619  CHKERR VecDestroy(&ghostFlux);
1621  }
1622 
1623  /**
1624  * \brief Calculate boundary conditions for fluxes
1625  * @return Error code
1626  */
1629  // clear vectors
1630  CHKERR VecZeroEntries(D0);
1631  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
1632  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
1633  // clear essential bc indices, it could have dofs from other mesh refinement
1634  bcIndices.clear();
1635  // set operator to calculate essential boundary conditions
1636  CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
1637  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
1638  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
1639  CHKERR VecAssemblyBegin(D0);
1640  CHKERR VecAssemblyEnd(D0);
1641  double norm2D0;
1642  CHKERR VecNorm(D0, NORM_2, &norm2D0);
1643  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1644  PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
1646  }
1647 
1648  /**
1649  * \brief Calculate inital pressure head distribution
1650  * @return Error code
1651  */
1654  // clear vectors
1655  CHKERR VecZeroEntries(D1);
1656  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
1657  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
1658  // Calculate initial pressure head on each element
1660  // Assemble vector
1661  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_REVERSE);
1662  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_REVERSE);
1663  CHKERR VecAssemblyBegin(D1);
1664  CHKERR VecAssemblyEnd(D1);
1665  // Calculate norm
1666  double norm2D1;
1667  CHKERR VecNorm(D1, NORM_2, &norm2D1);
1668  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1669  PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
1671  }
1672 
1673  /**
1674  * \brief solve problem
1675  * @return error code
1676  */
1677  MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
1679  if (set_initial_pc) {
1680  // Set initial head
1681  CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
1682  }
1683 
1684  // Initiate vector from data on the mesh
1685  CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
1686 
1687  // Create time solver
1688  TS ts;
1689  CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
1690  // Use backward Euler method
1691  CHKERR TSSetType(ts, TSBEULER);
1692  // Set final time
1693  double ftime = 1;
1694  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
1695  // Setup solver from commabd line
1696  CHKERR TSSetFromOptions(ts);
1697  // Set DM to TS
1698  CHKERR TSSetDM(ts, dM);
1699 #if PETSC_VERSION_GE(3, 7, 0)
1700  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1701 #endif
1702  // Set-up monitor
1703  TsCtx *ts_ctx;
1704  DMMoFEMGetTsCtx(dM, &ts_ctx);
1705  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1706 
1707  // This add SNES monitor, to show error by fields. It is dirty trick
1708  // to add monitor, so code is hidden from doxygen
1709  CHKERR TSSetSolution(ts, D);
1710  CHKERR TSSetUp(ts);
1711  SNES snes;
1712  CHKERR TSGetSNES(ts, &snes);
1713 
1714 #if PETSC_VERSION_GE(3, 7, 0)
1715  {
1716  PetscViewerAndFormat *vf;
1717  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1718  PETSC_VIEWER_DEFAULT, &vf);
1719  CHKERR SNESMonitorSet(
1720  snes,
1721  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1722  void *))SNESMonitorFields,
1723  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1724  }
1725 #else
1726  {
1727  CHKERR SNESMonitorSet(snes,
1728  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1729  void *))SNESMonitorFields,
1730  0, 0);
1731  }
1732 #endif
1733 
1734  CHKERR TSSolve(ts, D);
1735 
1736  // Get statistic form TS and print it
1737  CHKERR TSGetTime(ts, &ftime);
1738  PetscInt steps, snesfails, rejects, nonlinits, linits;
1739  CHKERR TSGetTimeStepNumber(ts, &steps);
1740  CHKERR TSGetSNESFailures(ts, &snesfails);
1741  CHKERR TSGetStepRejections(ts, &rejects);
1742  CHKERR TSGetSNESIterations(ts, &nonlinits);
1743  CHKERR TSGetKSPIterations(ts, &linits);
1744  PetscPrintf(PETSC_COMM_WORLD,
1745  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
1746  "%D, linits %D\n",
1747  steps, rejects, snesfails, ftime, nonlinits, linits);
1748 
1750  }
1751 };
1752 
1753 } // namespace MixTransport
1754 
1755 #endif // __UNSATURATD_FLOW_HPP__
static Index< 'K', 3 > K
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
EntitiesFieldData::EntData EntData
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
constexpr double alpha
@ ROW
Definition: definitions.h:136
@ MF_ZERO
Definition: definitions.h:111
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#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
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1061
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:758
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
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:481
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1080
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:811
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:229
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
double w(const double g, const double t)
double flux
impulse magnitude
double scale
Definition: plastic.cpp:103
Generic material model for unsaturated water transport.
virtual MoFEMErrorCode calC()
virtual MoFEMErrorCode calSe()
virtual void printMatParameters(const int id, const std::string &prefix) const =0
virtual MoFEMErrorCode calTheta()
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
static double ePsilon0
Regularization parameter.
virtual MoFEMErrorCode calK()
virtual MoFEMErrorCode calDiffC()
double sCale
Scale time dependent eq.
double K
Hydraulic conductivity [L/s].
Range tEts
Elements with this material.
double C
Capacity [S^2/L^2].
double Se
Effective saturation.
static double scaleZ
Scale z direction.
virtual double initialPcEval() const =0
Initialize head.
double h_t
rate of hydraulic head
static double ePsilon1
Regularization parameter.
virtual MoFEMErrorCode calDiffK()
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
VectorDouble valuesAtGaussPts
values at integration points on element
VectorDouble divergenceAtGaussPts
divergence at integration points on element
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
Class storing information about boundary condition.
boost::function< double(const double x, const double y, const double z)> hookFun
Set integration rule to boundary elements.
int operator()(int p_row, int p_col, int p_data) const
MoFEMErrorCode operator()()
function is run for every finite element
boost::shared_ptr< ForcesAndSourcesCore > fluxIntegrate
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
UnsaturatedFlowElement & cTx
MonitorPostProc(UnsaturatedFlowElement &ctx, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc, boost::shared_ptr< ForcesAndSourcesCore > flux_Integrate, const int frequency)
MoFEMErrorCode postProcess()
function is run at the end of loop
const int fRequency
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpDivTauU_HdivL2(UnsaturatedFlowElement &ctx, const std::string &flux_name_col, const std::string &val_name_row)
Constructor.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx, const std::string &val_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
OpIntegrateFluxes(UnsaturatedFlowElement &ctx, const std::string fluxes_name)
Constructor.
OpPostProcMaterial(UnsaturatedFlowElement &ctx, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
moab::Interface & postProcMesh
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
UnsaturatedFlowElement & cTx
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Evaluate boundary condition at the boundary.
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name, boost::shared_ptr< MethodForForceScaling > &value_scale)
Constructor.
boost::shared_ptr< MethodForForceScaling > valueScale
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrate boundary condition.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx, const std::string flux_name)
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx, const std::string &val_name_row, const std::string &flux_name_col)
Constructor.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Do calculations.
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Assemble matrix.
Set integration rule to volume elements.
int operator()(int, int, int p_data) const
Post proces method for volume element Assemble vectors and matrices and apply essential boundary cond...
postProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
boost::shared_ptr< ForcesAndSourcesCore > fePtr
Pre-peprocessing Set head pressute rate and get inital essential boundary conditions.
boost::shared_ptr< ForcesAndSourcesCore > fePtr
preProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
Implementation of operators, problem and finite elements for unsaturated flow.
MoFEMErrorCode calculateEssentialBc()
Calculate boundary conditions for fluxes.
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
Assemble tangent matrix.
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
Get value on boundary.
std::map< int, boost::shared_ptr< GenericMaterial > > MaterialsDoubleMap
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
Create finite element instances.
boost::shared_ptr< VectorDouble > headRateAtGaussPts
Vector keeps head rate.
boost::shared_ptr< MethodForForceScaling > scaleMethodFlux
Method scaling fluxes.
MoFEMErrorCode buildProblem(Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
Build problem.
Vec D1
Vector with inital head capillary pressure.
map< int, boost::shared_ptr< BcData > > BcMap
MaterialsDoubleMap dMatMap
materials database
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
Face element apply natural bc.
Vec ghostFlux
Ghost Vector of integrated fluxes.
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
add fields
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
Assemble residual vector.
UnsaturatedFlowElement(MoFEM::Interface &m_field)
boost::shared_ptr< FEMethod > tsMonitor
BcMap bcValueMap
Store boundary condition for head capillary pressure.
MoFEMErrorCode solveProblem(bool set_initial_pc=true)
solve problem
MoFEMErrorCode destroyMatrices()
Delete matrices and vector when no longer needed.
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.
MoFEMErrorCode calculateInitialPc()
Calculate inital pressure head distribution.
DM dM
Discrete manager for unsaturated flow problem.
MoFEMErrorCode createMatrices()
Create vectors and matrices.
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
Elemnet to calculate essential bc.
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
virtual MoFEMErrorCode getMaterial(const EntityHandle ent, int &block_id) const
For given element handle get material block Id.
boost::shared_ptr< MethodForForceScaling > scaleMethodValue
Method scaling values.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
structure for User Loop Methods on finite elements
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
TS ts
time solver
PetscReal ts_t
time
Vec & ts_F
residual vector
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Mat & ts_B
Preconditioner for ts_A.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Post processing.
Force scale operator for reading two columns.