v0.9.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  double sCale; ///< Scale time dependent eq.
42  // double Ks; ///< Saturated hydraulic conductivity [m/day]
43 
44  double h; ///< hydraulic head
45  double h_t; ///< rate of hydraulic head
46  // double h_flux_residual; ///< residual at point
47  double K; ///< Hydraulic conductivity [L/s]
48  double diffK; ///< Derivative of hydraulic conductivity [L/s * L^2/F]
49  double C; ///< Capacity [S^2/L^2]
50  double diffC; ///< Derivative of capacity [S^2/L^2 * L^2/F ]
51  double tHeta; ///< Water content
52  double Se; ///< Effective saturation
53 
54  Range tEts; ///< Elements with this material
55 
56  double x, y, z; ///< in meters (L)
57 
58  /**
59  * \brief Initialize head
60  * @return value of head
61  */
62  virtual double initalPcEval() const = 0;
63  virtual void printMatParameters(const int id,
64  const std::string &prefix) const = 0;
65 
66  virtual MoFEMErrorCode calK() {
68  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
69  "Not implemented how to calculate hydraulic conductivity");
71  }
72 
75  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
76  "Not implemented how to calculate derivative of hydraulic "
77  "conductivity");
79  }
80 
81  virtual MoFEMErrorCode calC() {
83  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
84  "Not implemented how to calculate capacity");
86  }
87 
90  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
91  "Not implemented how to calculate capacity");
93  }
94 
97  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
98  "Not implemented how to calculate capacity");
100  }
101 
102  virtual MoFEMErrorCode calSe() {
104  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
105  "Not implemented how to calculate capacity");
107  }
108 };
109 
110 /**
111  * \brief Implementation of operators, problem and finite elements for
112  * unsaturated flow
113  */
115 
116  DM dM; ///< Discrete manager for unsaturated flow problem
117 
119  : MixTransportElement(m_field), dM(PETSC_NULL), lastEvalBcValEnt(0),
121  lastEvalBcBlockFluxId(-1) {}
122 
124  if (dM != PETSC_NULL) {
125  CHKERR DMDestroy(&dM);
126  CHKERRABORT(PETSC_COMM_WORLD, ierr);
127  }
128  }
129 
130  typedef std::map<int, boost::shared_ptr<GenericMaterial>> MaterialsDoubleMap;
131  MaterialsDoubleMap dMatMap; ///< materials database
132 
133  /**
134  * \brief For given element handle get material block Id
135  * @param ent finite element entity handle
136  * @param block_id reference to returned block id
137  * @return error code
138  */
140  int &block_id) const {
142  for (MaterialsDoubleMap::const_iterator mit = dMatMap.begin();
143  mit != dMatMap.end(); mit++) {
144  if (mit->second->tEts.find(ent) != mit->second->tEts.end()) {
145  block_id = mit->first;
147  }
148  }
150  "Element not found, no material data");
152  }
153 
154  /**
155  * \brief Class storing information about boundary condition
156  */
157  struct BcData {
158  Range eNts;
159  double fixValue;
160  boost::function<double(const double x, const double y, const double z)>
162  BcData() : hookFun(NULL) {}
163  };
164  typedef map<int, boost::shared_ptr<BcData>> BcMap;
165  BcMap bcValueMap; ///< Store boundary condition for head capillary pressure
166 
169 
170  /**
171  * \brief Get value on boundary
172  * @param ent entity handle
173  * @param gg number of integration point
174  * @param x x-coordinate
175  * @param y y-coordinate
176  * @param z z-coordinate
177  * @param value returned value
178  * @return error code
179  */
180  MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg,
181  const double x, const double y, const double z,
182  double &value) {
184  int block_id = -1;
185  if (lastEvalBcValEnt == ent) {
186  block_id = lastEvalBcBlockValId;
187  } else {
188  for (BcMap::iterator it = bcValueMap.begin(); it != bcValueMap.end();
189  it++) {
190  if (it->second->eNts.find(ent) != it->second->eNts.end()) {
191  block_id = it->first;
192  }
193  }
194  lastEvalBcValEnt = ent;
195  lastEvalBcBlockValId = block_id;
196  }
197  if (block_id >= 0) {
198  if (bcValueMap.at(block_id)->hookFun) {
199  value = bcValueMap.at(block_id)->hookFun(x, y, z);
200  } else {
201  value = bcValueMap.at(block_id)->fixValue;
202  }
203  } else {
204  value = 0;
205  }
207  }
208 
212 
213  /**
214  * \brief essential (Neumann) boundary condition (set fluxes)
215  * @param ent handle to finite element entity
216  * @param x coord
217  * @param y coord
218  * @param z coord
219  * @param flux reference to flux which is set by function
220  * @return [description]
221  */
222  MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
223  const double y, const double z, double &flux) {
225  int block_id = -1;
226  if (lastEvalBcFluxEnt == ent) {
227  block_id = lastEvalBcBlockFluxId;
228  } else {
229  for (BcMap::iterator it = bcFluxMap.begin(); it != bcFluxMap.end();
230  it++) {
231  if (it->second->eNts.find(ent) != it->second->eNts.end()) {
232  block_id = it->first;
233  }
234  }
235  lastEvalBcFluxEnt = ent;
236  lastEvalBcBlockFluxId = block_id;
237  }
238  if (block_id >= 0) {
239  if (bcFluxMap.at(block_id)->hookFun) {
240  flux = bcFluxMap.at(block_id)->hookFun(x, y, z);
241  } else {
242  flux = bcFluxMap.at(block_id)->fixValue;
243  }
244  } else {
245  flux = 0;
246  }
248  }
249 
250  /**
251  * \brief Evaluate boundary condition at the boundary
252  */
255 
257  boost::shared_ptr<MethodForForceScaling> valueScale;
258 
259  /**
260  * \brief Constructor
261  */
262  OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name,
263  boost::shared_ptr<MethodForForceScaling> &value_scale)
265  fluxes_name, UserDataOperator::OPROW),
266  cTx(ctx), valueScale(value_scale) {}
267 
268  VectorDouble nF; ///< Vector of residuals
269 
270  /**
271  * \brief Integrate boundary condition
272  * @param side local index of entity
273  * @param type type of entity
274  * @param data data on entity
275  * @return error code
276  */
277  MoFEMErrorCode doWork(int side, EntityType type,
280  if (data.getFieldData().size() == 0)
282  // Get EntityHandle of the finite element
283  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
284  // Resize and clear vector
285  nF.resize(data.getIndices().size());
286  nF.clear();
287  // Get number of integration points
288  int nb_gauss_pts = data.getN().size1();
289  for (int gg = 0; gg < nb_gauss_pts; gg++) {
290  double x, y, z;
291  x = getCoordsAtGaussPts()(gg, 0);
292  y = getCoordsAtGaussPts()(gg, 1);
293  z = getCoordsAtGaussPts()(gg, 2);
294  double value;
295  // get value of boundary condition
296  CHKERR cTx.getBcOnValues(fe_ent, gg, x, y, z, value);
297  const double w = getGaussPts()(2, gg) * 0.5;
298  const double beta = w * (value - z);
299  noalias(nF) += beta * prod(data.getVectorN<3>(gg), getNormal());
300  }
301  // Scale vector if history evaluating method is given
302  Vec f = getFEMethod()->ts_F;
303  if (valueScale) {
304  CHKERR valueScale->scaleNf(getFEMethod(), nF);
305  }
306  // Assemble vector
307  CHKERR VecSetValues(f, data.getIndices().size(), &data.getIndices()[0],
308  &nF[0], ADD_VALUES);
310  }
311  };
312 
313  /**
314  * \brief Assemble flux residual
315  */
318 
320 
321  OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
323  flux_name, UserDataOperator::OPROW),
324  cTx(ctx) {}
325 
328 
329  MoFEMErrorCode doWork(int side, EntityType type,
332  const int nb_dofs = data.getIndices().size();
333  if (nb_dofs == 0)
335  nF.resize(nb_dofs, false);
336  nF.clear();
337  // Get EntityHandle of the finite element
338  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
339  // Get material block id
340  int block_id;
341  CHKERR cTx.getMaterial(fe_ent, block_id);
342  // Get material block
343  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
344  // block_data->printMatParameters(block_id,"Read material");
345  // Get base function
346  auto t_n_hdiv = data.getFTensor1N<3>();
347  // Get pressure
349  // Get flux
350  auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
351  // Coords at integration points
352  auto t_coords = getFTensor1CoordsAtGaussPts();
353  // Get integration weight
354  auto t_w = getFTensor0IntegrationWeight();
355  // Get volume
356  double vol = getVolume();
357  // Get material parameters
358  int nb_gauss_pts = data.getN().size1();
359  for (int gg = 0; gg != nb_gauss_pts; gg++) {
360  // Get divergence
361  CHKERR getDivergenceOfHDivBaseFunctions(side, type, data, gg, divVec);
362  const double alpha = t_w * vol;
363  block_data->h = t_h;
364  block_data->x = t_coords(0);
365  block_data->y = t_coords(1);
366  block_data->z = t_coords(2);
367  CHKERR block_data->calK();
368  const double K = block_data->K;
369  const double z = t_coords(2); /// z-coordinate at Gauss pt
370  // Calculate pressure gradient
371  noalias(nF) -= alpha * (t_h - z) * divVec;
372  // Calculate presure gradient from flux
373  FTensor::Tensor0<double *> t_nf(&*nF.begin());
374  for (int rr = 0; rr != nb_dofs; rr++) {
375  t_nf += alpha * (1 / K) * (t_n_hdiv(i) * t_flux(i));
376  ++t_n_hdiv; // move to next base function
377  ++t_nf; // move to next element in vector
378  }
379  ++t_h; // move to next integration point
380  ++t_flux;
381  ++t_coords;
382  ++t_w;
383  }
384  // Assemble residual
385  CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs,
386  &*data.getIndices().begin(), &*nF.begin(),
387  ADD_VALUES);
389  }
390  };
391 
392  /**
393  * Assemble mass residual
394  */
397 
399 
400  OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
402  val_name, UserDataOperator::OPROW),
403  cTx(ctx) {}
404 
406 
407  MoFEMErrorCode doWork(int side, EntityType type,
410  const int nb_dofs = data.getIndices().size();
411  if (nb_dofs == 0)
413  // Resize local element vector
414  nF.resize(nb_dofs, false);
415  nF.clear();
416  // Get EntityHandle of the finite element
417  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
418  // Get material block id
419  int block_id;
420  CHKERR cTx.getMaterial(fe_ent, block_id);
421  // Get material block
422  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
423  // Get pressure
425  // Get pressure rate
427  // Flux divergence
428  auto t_div_flux = getFTensor0FromVec(cTx.divergenceAtGaussPts);
429  // Get integration weight
430  auto t_w = getFTensor0IntegrationWeight();
431  // Coords at integration points
432  auto t_coords = getFTensor1CoordsAtGaussPts();
433  // Scale eq.
434  const double scale = block_data->sCale;
435  // Get volume
436  const double vol = getVolume();
437  // Get number of integration points
438  int nb_gauss_pts = data.getN().size1();
439  for (int gg = 0; gg != nb_gauss_pts; gg++) {
440  const double alpha = t_w * vol * scale;
441  block_data->h = t_h;
442  block_data->x = t_coords(0);
443  block_data->y = t_coords(1);
444  block_data->z = t_coords(2);
445  CHKERR block_data->calC();
446  const double C = block_data->C;
447  // Calculate flux conservation
448  noalias(nF) += (alpha * (t_div_flux + C * t_h_t)) * data.getN(gg);
449  ++t_h;
450  ++t_h_t;
451  ++t_div_flux;
452  ++t_coords;
453  ++t_w;
454  }
455  // Assemble local vector
456  Vec f = getFEMethod()->ts_F;
457  CHKERR VecSetValues(f, nb_dofs, &*data.getIndices().begin(), &*nF.begin(),
458  ADD_VALUES);
460  }
461  };
462 
465 
467 
469  const std::string flux_name)
471  flux_name, flux_name, UserDataOperator::OPROWCOL),
472  cTx(ctx) {
473  sYmm = true;
474  }
475 
477 
479 
480  /**
481  * \brief Assemble matrix
482  * @param row_side local index of row entity on element
483  * @param col_side local index of col entity on element
484  * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
485  * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
486  * @param row_data data for row
487  * @param col_data data for col
488  * @return error code
489  */
490  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
491  EntityType col_type,
495  const int nb_row = row_data.getIndices().size();
496  const int nb_col = col_data.getIndices().size();
497  if (nb_row == 0)
499  if (nb_col == 0)
501  nN.resize(nb_row, nb_col, false);
502  nN.clear();
503  // Get EntityHandle of the finite element
504  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
505  // Get material block id
506  int block_id;
507  CHKERR cTx.getMaterial(fe_ent, block_id);
508  // Get material block
509  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
510  // Get pressure
512  // Coords at integration points
513  auto t_coords = getFTensor1CoordsAtGaussPts();
514  // Get base functions
515  auto t_n_hdiv_row = row_data.getFTensor1N<3>();
516  // Get integration weight
517  auto t_w = getFTensor0IntegrationWeight();
518  // Get volume
519  const double vol = getVolume();
520  int nb_gauss_pts = row_data.getN().size1();
521  for (int gg = 0; gg != nb_gauss_pts; gg++) {
522  block_data->h = t_h;
523  block_data->x = t_coords(0);
524  block_data->y = t_coords(1);
525  block_data->z = t_coords(2);
526  CHKERR block_data->calK();
527  const double K = block_data->K;
528  // get integration weight and multiply by element volume
529  const double alpha = t_w * vol;
530  const double beta = alpha * (1 / K);
531  FTensor::Tensor0<double *> t_a(&*nN.data().begin());
532  for (int kk = 0; kk != nb_row; kk++) {
533  auto t_n_hdiv_col = col_data.getFTensor1N<3>(gg, 0);
534  for (int ll = 0; ll != nb_col; ll++) {
535  t_a += beta * (t_n_hdiv_row(j) * t_n_hdiv_col(j));
536  ++t_n_hdiv_col;
537  ++t_a;
538  }
539  ++t_n_hdiv_row;
540  }
541  ++t_coords;
542  ++t_h;
543  ++t_w;
544  }
545  Mat a = getFEMethod()->ts_B;
546  CHKERR MatSetValues(a, nb_row, &*row_data.getIndices().begin(), nb_col,
547  &*col_data.getIndices().begin(), &*nN.data().begin(),
548  ADD_VALUES);
549  // matrix is symmetric, assemble other part
550  if (row_side != col_side || row_type != col_type) {
551  transNN.resize(nb_col, nb_row);
552  noalias(transNN) = trans(nN);
553  CHKERR MatSetValues(a, nb_col, &*col_data.getIndices().begin(), nb_row,
554  &*row_data.getIndices().begin(),
555  &*transNN.data().begin(), ADD_VALUES);
556  }
558  }
559  };
560 
561  struct OpVU_L2L2
563 
565 
566  OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
568  value_name, value_name, UserDataOperator::OPROWCOL),
569  cTx(ctx) {
570  sYmm = true;
571  }
572 
574 
575  /**
576  * \brief Assemble matrix
577  * @param row_side local index of row entity on element
578  * @param col_side local index of col entity on element
579  * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
580  * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
581  * @param row_data data for row
582  * @param col_data data for col
583  * @return error code
584  */
585  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
586  EntityType col_type,
590  int nb_row = row_data.getIndices().size();
591  int nb_col = col_data.getIndices().size();
592  if (nb_row == 0)
594  if (nb_col == 0)
596  nN.resize(nb_row, nb_col, false);
597  nN.clear();
598  // Get EntityHandle of the finite element
599  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
600  // Get material block id
601  int block_id;
602  CHKERR cTx.getMaterial(fe_ent, block_id);
603  // Get material block
604  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
605  // Get pressure
607  // // Get pressure
608  // auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
609  // Get pressure rate
611  // Get integration weight
612  auto t_w = getFTensor0IntegrationWeight();
613  // Coords at integration points
614  auto t_coords = getFTensor1CoordsAtGaussPts();
615  // Scale eq.
616  const double scale = block_data->sCale;
617  // Time step factor
618  double ts_a = getFEMethod()->ts_a;
619  // get volume
620  const double vol = getVolume();
621  int nb_gauss_pts = row_data.getN().size1();
622  // get base functions on rows
623  auto t_n_row = row_data.getFTensor0N();
624  for (int gg = 0; gg != nb_gauss_pts; gg++) {
625  // get integration weight and multiply by element volume
626  double alpha = t_w * vol * scale;
627  // evaluate material model at integration points
628  // to calculate capacity and tangent of capacity term
629  block_data->h = t_h;
630  block_data->h_t = t_h_t;
631  block_data->x = t_coords(0);
632  block_data->y = t_coords(1);
633  block_data->z = t_coords(2);
634  CHKERR block_data->calC();
635  CHKERR block_data->calDiffC();
636  const double C = block_data->C;
637  const double diffC = block_data->diffC;
638  // assemble local entity tangent matrix
640  &*nN.data().begin());
641  // iterate base functions on rows
642  for (int kk = 0; kk != nb_row; kk++) {
643  // get first base function on column at integration point gg
644  auto t_n_col = col_data.getFTensor0N(gg, 0);
645  // iterate base functions on columns
646  for (int ll = 0; ll != nb_col; ll++) {
647  // assemble elements of local matrix
648  t_a += (alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
649  ++t_n_col; // move to next base function on column
650  ++t_a; // move to next element in local tangent matrix
651  }
652  ++t_n_row; // move to next base function on row
653  }
654  ++t_w; // move to next integration weight
655  ++t_coords; // move to next coordinate at integration point
656  ++t_h; // move to next capillary head at integration point
657  // ++t_flux_residual;
658  ++t_h_t; // move to next capillary head rate at integration point
659  }
660  Mat a = getFEMethod()->ts_B;
661  CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
662  &col_data.getIndices()[0], &*nN.data().begin(),
663  ADD_VALUES);
665  }
666  };
667 
670 
672 
673  /**
674  * \brief Constructor
675  */
677  const std::string &val_name_row,
678  const std::string &flux_name_col)
680  val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
681  cTx(ctx) {}
682 
685 
686  /**
687  * \brief Do calculations
688  * @param row_side local index of entity on row
689  * @param col_side local index of entity on column
690  * @param row_type type of row entity
691  * @param col_type type of col entity
692  * @param row_data row data structure carrying information about base
693  * functions, DOFs indices, etc.
694  * @param col_data column data structure carrying information about base
695  * functions, DOFs indices, etc.
696  * @return error code
697  */
698  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
699  EntityType col_type,
703  int nb_row = row_data.getFieldData().size();
704  int nb_col = col_data.getFieldData().size();
705  if (nb_row == 0)
707  if (nb_col == 0)
709  // Get EntityHandle of the finite element
710  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
711  // Get material block id
712  int block_id;
713  CHKERR cTx.getMaterial(fe_ent, block_id);
714  // Get material block
715  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
716  nN.resize(nb_row, nb_col, false);
717  divVec.resize(nb_col, false);
718  nN.clear();
719  // Scale eq.
720  const double scale = block_data->sCale;
721  int nb_gauss_pts = row_data.getN().size1();
722  for (int gg = 0; gg < nb_gauss_pts; gg++) {
723  double alpha = getGaussPts()(3, gg) * getVolume() * scale;
724  CHKERR getDivergenceOfHDivBaseFunctions(col_side, col_type, col_data,
725  gg, divVec);
726  noalias(nN) += alpha * outer_prod(row_data.getN(gg), divVec);
727  }
728  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
729  &row_data.getIndices()[0], nb_col,
730  &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
732  }
733  };
734 
737 
739 
740  /**
741  * \brief Constructor
742  */
744  const std::string &flux_name_col,
745  const std::string &val_name_row)
747  flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
748  cTx(ctx) {}
749 
753 
754  /**
755  * \brief Do calculations
756  * @param row_side local index of entity on row
757  * @param col_side local index of entity on column
758  * @param row_type type of row entity
759  * @param col_type type of col entity
760  * @param row_data row data structure carrying information about base
761  * functions, DOFs indices, etc.
762  * @param col_data column data structure carrying information about base
763  * functions, DOFs indices, etc.
764  * @return error code
765  */
766  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
767  EntityType col_type,
771  int nb_row = row_data.getFieldData().size();
772  int nb_col = col_data.getFieldData().size();
773  if (nb_row == 0)
775  if (nb_col == 0)
777  nN.resize(nb_row, nb_col, false);
778  divVec.resize(nb_row, false);
779  nN.clear();
780  // Get EntityHandle of the finite element
781  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
782  // Get material block id
783  int block_id;
784  CHKERR cTx.getMaterial(fe_ent, block_id);
785  // Get material block
786  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
787  // Get pressure
789  // Get flux
790  auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
791  // Coords at integration points
792  auto t_coords = getFTensor1CoordsAtGaussPts();
793  // Get integration weight
794  auto t_w = getFTensor0IntegrationWeight();
795  // Get base function
796  auto t_n_hdiv_row = row_data.getFTensor1N<3>();
797  // Get volume
798  double vol = getVolume();
799  int nb_gauss_pts = row_data.getN().size1();
800  for (int gg = 0; gg != nb_gauss_pts; gg++) {
801  block_data->h = t_h;
802  block_data->x = t_coords(0);
803  block_data->y = t_coords(1);
804  block_data->z = t_coords(2);
805  CHKERR block_data->calK();
806  CHKERR block_data->calDiffK();
807  const double K = block_data->K;
808  // const double z = t_coords(2);
809  const double KK = K * K;
810  const double diffK = block_data->diffK;
811  double alpha = t_w * vol;
812  CHKERR getDivergenceOfHDivBaseFunctions(row_side, row_type, row_data,
813  gg, divVec);
814  noalias(nN) -= alpha * outer_prod(divVec, col_data.getN(gg));
815  FTensor::Tensor0<double *> t_a(&*nN.data().begin());
816  for (int rr = 0; rr != nb_row; rr++) {
817  double beta = alpha * (-diffK / KK) * (t_n_hdiv_row(i) * t_flux(i));
818  auto t_n_col = col_data.getFTensor0N(gg, 0);
819  for (int cc = 0; cc != nb_col; cc++) {
820  t_a += beta * t_n_col;
821  ++t_n_col;
822  ++t_a;
823  }
824  ++t_n_hdiv_row;
825  }
826  ++t_w;
827  ++t_coords;
828  ++t_h;
829  ++t_flux;
830  }
831  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
832  &row_data.getIndices()[0], nb_col,
833  &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
835  }
836  };
837 
842  const std::string &val_name)
844  val_name, UserDataOperator::OPROW),
845  cTx(ctx) {}
846 
849 
850  MoFEMErrorCode doWork(int side, EntityType type,
853  if (data.getFieldData().size() == 0)
855  int nb_dofs = data.getFieldData().size();
856  int nb_gauss_pts = data.getN().size1();
857  if (nb_dofs != static_cast<int>(data.getN().size2())) {
858  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
859  "wrong number of dofs");
860  }
861  nN.resize(nb_dofs, nb_dofs, false);
862  nF.resize(nb_dofs, false);
863  nN.clear();
864  nF.clear();
865 
866  // Get EntityHandle of the finite element
867  EntityHandle fe_ent = getFEEntityHandle();
868  // Get material block id
869  int block_id;
870  CHKERR cTx.getMaterial(fe_ent, block_id);
871  // Get material block
872  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
873 
874  // loop over integration points
875  for (int gg = 0; gg < nb_gauss_pts; gg++) {
876  // get coordinates at integration point
877  block_data->x = getCoordsAtGaussPts()(gg, 0);
878  block_data->y = getCoordsAtGaussPts()(gg, 1);
879  block_data->z = getCoordsAtGaussPts()(gg, 2);
880  // get weight for integration rule
881  double alpha = getGaussPts()(2, gg) * getVolume();
882  nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
883  nF += alpha * block_data->initalPcEval() * data.getN(gg);
884  }
885 
886  // factor matrix
888  // solve local problem
889  cholesky_solve(nN, nF, ublas::lower());
890 
891  // set solution to vector
892  CHKERR VecSetValues(cTx.D1, nb_dofs, &*data.getIndices().begin(),
893  &*nF.begin(), INSERT_VALUES);
894 
896  }
897  };
898 
902 
903  /**
904  * \brief Constructor
905  */
907  const std::string fluxes_name)
909  fluxes_name, UserDataOperator::OPROW),
910  cTx(ctx) {}
911 
913 
914  /**
915  * \brief Integrate boundary condition
916  * @param side local index of entity
917  * @param type type of entity
918  * @param data data on entity
919  * @return error code
920  */
921  MoFEMErrorCode doWork(int side, EntityType type,
924  int nb_dofs = data.getFieldData().size();
925  if (nb_dofs == 0)
927  // Get base function
928  auto t_n_hdiv = data.getFTensor1N<3>();
929  // get normal of face
930  auto t_normal = getFTensor1Normal();
931  // Integration weight
932  auto t_w = getFTensor0IntegrationWeight();
933  double flux_on_entity = 0;
934  int nb_gauss_pts = data.getN().size1();
935  for (int gg = 0; gg < nb_gauss_pts; gg++) {
936  auto t_data = data.getFTensor0FieldData();
937  for (int rr = 0; rr != nb_dofs; rr++) {
938  flux_on_entity -= (0.5 * t_data * t_w) * (t_n_hdiv(i) * t_normal(i));
939  ++t_n_hdiv;
940  ++t_data;
941  }
942  ++t_w;
943  }
944  CHKERR VecSetValue(cTx.ghostFlux, 0, flux_on_entity, ADD_VALUES);
946  }
947  };
948 
949  /**
950  * Operator used to post-process results for unsaturated infiltration problem.
951  * Operator should with element for post-processing results, i.e.
952  * PostProcVolumeOnRefinedMesh
953  */
956 
958  moab::Interface &postProcMesh;
959  std::vector<EntityHandle> &mapGaussPts;
960 
962  moab::Interface &post_proc_mesh,
963  std::vector<EntityHandle> &map_gauss_pts,
964  const std::string field_name)
967  cTx(ctx), postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
968 
969  MoFEMErrorCode doWork(int side, EntityType type,
972  int nb_dofs = data.getFieldData().size();
973  if (nb_dofs == 0)
975 
976  // Get EntityHandle of the finite element
977  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
978  // Get material block id
979  int block_id;
980  CHKERR cTx.getMaterial(fe_ent, block_id);
981  // Get material block
982  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
983 
984  // Set bloc Id
985  Tag th_id;
986  int def_block_id = -1;
987  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
988  MB_TAG_CREAT | MB_TAG_SPARSE,
989  &def_block_id);
990  CHKERR postProcMesh.tag_set_data(th_id, &fe_ent, 1, &block_id);
991 
992  // Create mesh tag. Tags are created on post-processing mesh and
993  // visable in post-processor, e.g. Paraview
994  double zero = 0;
995  Tag th_theta;
996  CHKERR postProcMesh.tag_get_handle("THETA", 1, MB_TYPE_DOUBLE, th_theta,
997  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
998  Tag th_se;
999  CHKERR postProcMesh.tag_get_handle("Se", 1, MB_TYPE_DOUBLE, th_se,
1000  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1001  // Tag th_ks;
1002  // CHKERR postProcMesh.tag_get_handle(
1003  // "Ks",1,MB_TYPE_DOUBLE,th_ks,
1004  // MB_TAG_CREAT|MB_TAG_SPARSE,&zero
1005  // );
1006  // CHKERR postProcMesh.tag_set_data(th_ks,&fe_ent,1,&block_data->Ks);
1007  Tag th_k;
1008  CHKERR postProcMesh.tag_get_handle("K", 1, MB_TYPE_DOUBLE, th_k,
1009  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1010  Tag th_c;
1011  CHKERR postProcMesh.tag_get_handle("C", 1, MB_TYPE_DOUBLE, th_c,
1012  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1013 
1014  // Get pressure at integration points
1016  // Coords at integration points
1017  auto t_coords = getFTensor1CoordsAtGaussPts();
1018 
1019  int nb_gauss_pts = data.getN().size1();
1020  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1021  block_data->h = t_h;
1022  block_data->x = t_coords(0);
1023  block_data->y = t_coords(1);
1024  block_data->z = t_coords(2);
1025  // Calculate theta (water content) and save it on mesh tags
1026  CHKERR block_data->calTheta();
1027  double theta = block_data->tHeta;
1028  CHKERR postProcMesh.tag_set_data(th_theta, &mapGaussPts[gg], 1, &theta);
1029  CHKERR block_data->calSe();
1030  // Calculate Se (effective saturation and save it on the mesh tags)
1031  double Se = block_data->Se;
1032  CHKERR postProcMesh.tag_set_data(th_se, &mapGaussPts[gg], 1, &Se);
1033  // Calculate K (hydraulic conductivity) and save it on the mesh tags
1034  CHKERR block_data->calK();
1035  double K = block_data->K;
1036  CHKERR postProcMesh.tag_set_data(th_k, &mapGaussPts[gg], 1, &K);
1037  // Calculate water capacity and save it on the mesh tags
1038  CHKERR block_data->calC();
1039  double C = block_data->C;
1040  CHKERR postProcMesh.tag_set_data(th_c, &mapGaussPts[gg], 1, &C);
1041  ++t_h;
1042  ++t_coords;
1043  }
1044 
1046  }
1047  };
1048 
1049  /**
1050  * Finite element implementation called by TS monitor. Element calls other
1051  * finite elements to evaluate material properties and save results on the
1052  * mesh.
1053  *
1054  * \note Element overloaded only FEMethod::postProcess methos where other
1055  * elements are called.
1056  */
1058 
1060  boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProc;
1061  boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
1062 
1063  const int fRequency;
1064 
1066  boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc,
1067  boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
1068  const int frequency)
1069  : cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
1070  fRequency(frequency) {}
1071 
1075  }
1076 
1080  }
1081 
1084 
1085  // Get time step
1086  int step;
1087  CHKERR TSGetTimeStepNumber(ts, &step);
1088 
1089  if ((step) % fRequency == 0) {
1090  // Post-process results and save in the file
1091  PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
1092  fRequency);
1094  CHKERR postProc->writeFile(
1095  string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
1096  }
1097 
1098  // Integrate fluxes on faces where pressure head is applied
1099  CHKERR VecZeroEntries(cTx.ghostFlux);
1100  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1101  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1102  // Run finite element to integrate fluxes
1104  CHKERR VecAssemblyBegin(cTx.ghostFlux);
1105  CHKERR VecAssemblyEnd(cTx.ghostFlux);
1106  // accumulate errors from processors
1107  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1108  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1109  // scatter errors to all processors
1110  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1111  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1112  double *ghost_flux;
1113  CHKERR VecGetArray(cTx.ghostFlux, &ghost_flux);
1114  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Flux at time %6.4g %6.4g\n", ts_t,
1115  ghost_flux[0]);
1116  CHKERR VecRestoreArray(cTx.ghostFlux, &ghost_flux);
1117 
1119  }
1120  };
1121 
1122  /// \brief add fields
1123  MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
1124  const int order) {
1126  // Fields
1129  CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
1130  // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
1131 
1132  // meshset consisting all entities in mesh
1133  EntityHandle root_set = mField.get_moab().get_root_set();
1134  // add entities to field
1135 
1137  if (it->getName().compare(0, 4, "SOIL") != 0)
1138  continue;
1139  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1140  MBTET, fluxes);
1141  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1142  MBTET, values);
1143  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1144  MBTET, values + "_t");
1145  // CHKERR mField.add_ents_to_field_by_type(
1146  // dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
1147  // );
1148  }
1149 
1150  CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
1151  CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
1152  CHKERR mField.set_field_order(root_set, MBTET, values, order);
1153  CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
1154  // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
1156  }
1157 
1158  /// \brief add finite elements
1159  MoFEMErrorCode addFiniteElements(const std::string &fluxes_name,
1160  const std::string &values_name) {
1162 
1163  // Define element "MIX". Note that this element will work with fluxes_name
1164  // and values_name. This reflect bilinear form for the problem
1173  values_name + "_t");
1174  // CHKERR
1175  // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
1176 
1178  if (it->getName().compare(0, 4, "SOIL") != 0)
1179  continue;
1181  dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
1182  }
1183 
1184  // Define element to integrate natural boundary conditions, i.e. set values.
1185  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
1187  fluxes_name);
1189  fluxes_name);
1191  fluxes_name);
1193  values_name);
1194 
1196  if (it->getName().compare(0, 4, "HEAD") != 0)
1197  continue;
1199  bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
1200  }
1201 
1202  // Define element to apply essential boundary conditions.
1203  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
1205  fluxes_name);
1207  fluxes_name);
1209  fluxes_name);
1211  values_name);
1212 
1214  if (it->getName().compare(0, 4, "FLUX") != 0)
1215  continue;
1217  bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
1218  }
1219 
1221  }
1222 
1223  /**
1224  * \brief Build problem
1225  * @param ref_level mesh refinement on which mesh problem you like to built.
1226  * @return error code
1227  */
1230 
1231  // Build fields
1233  // Build finite elements
1235  CHKERR mField.build_finite_elements("MIX_BCFLUX");
1236  CHKERR mField.build_finite_elements("MIX_BCVALUE");
1237  // Build adjacencies of degrees of freedom and elements
1238  CHKERR mField.build_adjacencies(ref_level);
1239 
1240  // create DM instance
1241  CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
1242  // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
1243  CHKERR DMSetType(dM, "DMMOFEM");
1244  // mesh is portioned, each process keeps only part of problem
1245  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1246  // creates problem in DM
1247  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
1248  // discretised problem creates square matrix (that makes some optimizations)
1249  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1250  // set DM options from command line
1251  CHKERR DMSetFromOptions(dM);
1252  // add finite elements
1253  CHKERR DMMoFEMAddElement(dM, "MIX");
1254  CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
1255  CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
1256  // constructor data structures
1257  CHKERR DMSetUp(dM);
1258 
1259  PetscSection section;
1260  CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
1261  CHKERR DMSetDefaultSection(dM, section);
1262  CHKERR DMSetDefaultGlobalSection(dM, section);
1263  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
1264  CHKERR PetscSectionDestroy(&section);
1265 
1267  }
1268 
1269  boost::shared_ptr<ForcesAndSourcesCore>
1270  feFaceBc; ///< Elemnet to calculate essential bc
1271  boost::shared_ptr<ForcesAndSourcesCore>
1272  feFaceRhs; ///< Face element apply natural bc
1273  boost::shared_ptr<ForcesAndSourcesCore>
1274  feVolInitialPc; ///< Calculate inital boundary conditions
1275  boost::shared_ptr<ForcesAndSourcesCore>
1276  feVolRhs; ///< Assemble residual vector
1277  boost::shared_ptr<ForcesAndSourcesCore> feVolLhs; ///< Assemble tangent matrix
1278  boost::shared_ptr<MethodForForceScaling>
1279  scaleMethodFlux; ///< Method scaling fluxes
1280  boost::shared_ptr<MethodForForceScaling>
1281  scaleMethodValue; ///< Method scaling values
1282  boost::shared_ptr<FEMethod> tsMonitor; ///< Element used by TS monitor to
1283  ///< postprocess results at time step
1284 
1285  boost::shared_ptr<VectorDouble>
1286  headRateAtGaussPts; ///< Vector keeps head rate
1287  // boost::shared_ptr<VectorDouble> resAtGaussPts; ///< Residual field
1288 
1289  /**
1290  * \brief Set integration rule to volume elements
1291  *
1292  */
1293  struct VolRule {
1294  int operator()(int, int, int p_data) const { return 2 * p_data + p_data; }
1295  };
1296  /**
1297  * \brief Set integration rule to boundary elements
1298  *
1299  */
1300  struct FaceRule {
1301  int operator()(int p_row, int p_col, int p_data) const {
1302  return 2 * p_data;
1303  }
1304  };
1305 
1306  std::vector<int> bcVecIds;
1308 
1309  /**
1310  * \brief Pre-peprocessing
1311  * Set head pressute rate and get inital essential boundary conditions
1312  */
1313  struct preProcessVol {
1315  boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1316  // std::string mArk;
1317 
1320  boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
1321  )
1322  : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1325  // Update pressure rates
1326  CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1327  fePtr->problemPtr, "VALUES", string("VALUES") + "_t", ROW,
1328  fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1329  switch (fePtr->ts_ctx) {
1330  case TSMethod::CTX_TSSETIFUNCTION:
1331  if (!cTx.bcIndices.empty()) {
1332  double scale;
1333  CHKERR cTx.scaleMethodFlux->getForceScale(fePtr->ts_t, scale);
1334  if (cTx.bcVecIds.size() != cTx.bcIndices.size()) {
1335  cTx.bcVecIds.insert(cTx.bcVecIds.begin(), cTx.bcIndices.begin(),
1336  cTx.bcIndices.end());
1337  cTx.bcVecVals.resize(cTx.bcVecIds.size(), false);
1338  cTx.vecValsOnBc.resize(cTx.bcVecIds.size(), false);
1339  }
1340  CHKERR VecGetValues(cTx.D0, cTx.bcVecIds.size(),
1341  &*cTx.bcVecIds.begin(), &*cTx.bcVecVals.begin());
1342  CHKERR VecGetValues(fePtr->ts_u, cTx.bcVecIds.size(),
1343  &*cTx.bcVecIds.begin(),
1344  &*cTx.vecValsOnBc.begin());
1345  cTx.bcVecVals *= scale;
1346  VectorDouble::iterator vit = cTx.bcVecVals.begin();
1347  const NumeredDofEntity *dof_ptr;
1348  for (std::vector<int>::iterator it = cTx.bcVecIds.begin();
1349  it != cTx.bcVecIds.end(); it++, vit++) {
1350  CHKERR fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it,
1351  &dof_ptr);
1352  dof_ptr->getFieldData() = *vit;
1353  }
1354  } else {
1355  cTx.bcVecIds.resize(0);
1356  cTx.bcVecVals.resize(0);
1357  cTx.vecValsOnBc.resize(0);
1358  }
1359  break;
1360  default:
1361  // don nothing
1362  break;
1363  }
1365  }
1366  };
1367 
1368  /**
1369  * \brief Post proces method for volume element
1370  * Assemble vectors and matrices and apply essential boundary conditions
1371  */
1374  boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1375  // std::string mArk;
1378  boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
1379  )
1380  : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1383  switch (fePtr->ts_ctx) {
1384  case TSMethod::CTX_TSSETIJACOBIAN: {
1385  CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1386  CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1387  // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1388  // std::string wait;
1389  // std::cin >> wait;
1390  CHKERR MatZeroRowsColumns(fePtr->ts_B, cTx.bcVecIds.size(),
1391  &*cTx.bcVecIds.begin(), 1, PETSC_NULL,
1392  PETSC_NULL);
1393  CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1394  CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1395  // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1396  // std::string wait;
1397  // std::cin >> wait;
1398  } break;
1399  case TSMethod::CTX_TSSETIFUNCTION: {
1400  CHKERR VecAssemblyBegin(fePtr->ts_F);
1401  CHKERR VecAssemblyEnd(fePtr->ts_F);
1402  if (!cTx.bcVecIds.empty()) {
1404  // cerr << mArk << endl;
1405  // cerr << "a " << cTx.vecValsOnBc << endl;
1406  // cerr << "a " << cTx.bcVecVals << endl;
1407  CHKERR VecSetValues(fePtr->ts_F, cTx.bcVecIds.size(),
1408  &*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
1409  INSERT_VALUES);
1410  }
1411  CHKERR VecAssemblyBegin(fePtr->ts_F);
1412  CHKERR VecAssemblyEnd(fePtr->ts_F);
1413  // CHKERR VecView(fePtr->ts_F,PETSC_VIEWER_STDOUT_WORLD);
1414  // CHKERR
1415  // fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1416  // fePtr->problemPtr,"VALUES",string("FLUXES")+"_residual",
1417  // ROW,fePtr->ts_F,INSERT_VALUES,SCATTER_REVERSE
1418  // );
1419  } break;
1420  default:
1421  // don nothing
1422  break;
1423  }
1425  }
1426  };
1427 
1428  /**
1429  * \brief Create finite element instances
1430  * @param vol_rule integration rule for volume element
1431  * @param face_rule integration rule for boundary element
1432  * @return error code
1433  */
1435  setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule = VolRule(),
1436  ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule()) {
1438 
1439  // create finite element instances
1440  feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
1442  feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1444  feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
1446  feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
1448  feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1450  // set integration rule to elements
1451  feFaceBc->getRuleHook = face_rule;
1452  feFaceRhs->getRuleHook = face_rule;
1453  feVolInitialPc->getRuleHook = vol_rule;
1454  feVolLhs->getRuleHook = vol_rule;
1455  feVolRhs->getRuleHook = vol_rule;
1456  // set function hook for finite element postprocessing stage
1457  feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
1458  feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
1459  feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
1460  feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
1461 
1462  // create method for setting history for fluxes on boundary
1463  scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
1464  new TimeForceScale("-flux_history", false));
1465 
1466  // create method for setting history for presure heads on boundary
1467  scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
1468  new TimeForceScale("-head_history", false));
1469 
1470  // Set operator to calculate essential boundary conditions
1471  feFaceBc->getOpPtrVector().push_back(
1472  new OpEvaluateBcOnFluxes(*this, "FLUXES"));
1473 
1474  // Set operator to calculate initial capillary pressure
1475  feVolInitialPc->getOpPtrVector().push_back(
1476  new OpEvaluateInitiallHead(*this, "VALUES"));
1477 
1478  // set residual face from Neumann terms, i.e. applied pressure
1479  feFaceRhs->getOpPtrVector().push_back(
1480  new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
1481  // set residual finite element operators
1482  headRateAtGaussPts = boost::make_shared<VectorDouble>();
1483  // resAtGaussPts = boost::make_shared<VectorDouble>();
1484  feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1485  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1486  feVolRhs->getOpPtrVector().push_back(
1487  new OpValuesAtGaussPts(*this, "VALUES"));
1488  feVolRhs->getOpPtrVector().push_back(
1489  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1490  feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
1491  feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
1492  feVolRhs->getOpPtrVector().back().opType =
1493  ForcesAndSourcesCore::UserDataOperator::OPROW;
1494  // set tangent matrix finite element operators
1495  feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1496  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1497  // feVolLhs->getOpPtrVector().push_back(
1498  // new
1499  // OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
1500  // );
1501  feVolLhs->getOpPtrVector().push_back(
1502  new OpValuesAtGaussPts(*this, "VALUES"));
1503  feVolLhs->getOpPtrVector().push_back(
1504  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1505  feVolLhs->getOpPtrVector().push_back(
1506  new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
1507  feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
1508  feVolLhs->getOpPtrVector().push_back(
1509  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
1510  feVolLhs->getOpPtrVector().push_back(
1511  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
1512 
1513  // Adding finite elements to DM, time solver will ask for it to assemble
1514  // tangent matrices and residuals
1515  boost::shared_ptr<FEMethod> null;
1516  CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
1517  CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
1518  CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
1519 
1520  // setting up post-processing
1521  boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_process(
1523  CHKERR post_process->generateReferenceElementMesh();
1524  CHKERR post_process->addFieldValuesPostProc("VALUES");
1525  CHKERR post_process->addFieldValuesPostProc("VALUES_t");
1526  // CHKERR post_process->addFieldValuesPostProc("FLUXES_residual");
1527  CHKERR post_process->addFieldValuesPostProc("FLUXES");
1528  post_process->getOpPtrVector().push_back(
1529  new OpValuesAtGaussPts(*this, "VALUES"));
1530  post_process->getOpPtrVector().push_back(
1531  new OpPostProcMaterial(*this, post_process->postProcMesh,
1532  post_process->mapGaussPts, "VALUES"));
1533 
1534  // Integrate fluxes on boundary
1535  boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
1536  flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
1538  flux_integrate->getOpPtrVector().push_back(
1539  new OpIntegrateFluxes(*this, "FLUXES"));
1540  int frequency = 1;
1541  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
1542  CHKERR PetscOptionsInt(
1543  "-how_often_output",
1544  "frequency how often results are dumped on hard disk", "", frequency,
1545  &frequency, NULL);
1546  ierr = PetscOptionsEnd();
1547  CHKERRG(ierr);
1548 
1549  tsMonitor = boost::shared_ptr<FEMethod>(
1550  new MonitorPostProc(*this, post_process, flux_integrate, frequency));
1551  TsCtx *ts_ctx;
1552  DMMoFEMGetTsCtx(dM, &ts_ctx);
1553  ts_ctx->get_postProcess_to_do_Monitor().push_back(tsMonitor);
1555  }
1556 
1557  Vec D1; ///< Vector with inital head capillary pressure
1558  Vec ghostFlux; ///< Ghost Vector of integrated fluxes
1559 
1560  /**
1561  * \brief Create vectors and matrices
1562  * @return Error code
1563  */
1566  CHKERR DMCreateMatrix(dM, &Aij);
1567  CHKERR DMCreateGlobalVector(dM, &D0);
1568  CHKERR VecDuplicate(D0, &D1);
1569  CHKERR VecDuplicate(D0, &D);
1570  CHKERR VecDuplicate(D0, &F);
1571  int ghosts[] = {0};
1572  int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
1573  int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
1574  CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
1575  &ghostFlux);
1577  }
1578 
1579  /**
1580  * \brief Delete matrices and vector when no longer needed
1581  * @return error code
1582  */
1585  CHKERR MatDestroy(&Aij);
1586  CHKERR VecDestroy(&D);
1587  CHKERR VecDestroy(&D0);
1588  CHKERR VecDestroy(&D1);
1589  CHKERR VecDestroy(&F);
1590  CHKERR VecDestroy(&ghostFlux);
1592  }
1593 
1594  /**
1595  * \brief Calculate boundary conditions for fluxes
1596  * @return Error code
1597  */
1600  // clear vectors
1601  CHKERR VecZeroEntries(D0);
1602  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
1603  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
1604  // clear essential bc indices, it could have dofs from other mesh refinement
1605  bcIndices.clear();
1606  // set operator to calculate essential boundary conditions
1607  CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
1608  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
1609  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
1610  CHKERR VecAssemblyBegin(D0);
1611  CHKERR VecAssemblyEnd(D0);
1612  double norm2D0;
1613  CHKERR VecNorm(D0, NORM_2, &norm2D0);
1614  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1615  PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
1617  }
1618 
1619  /**
1620  * \brief Calculate inital pressure head distribution
1621  * @return Error code
1622  */
1625  // clear vectors
1626  CHKERR VecZeroEntries(D1);
1627  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
1628  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
1629  // Calculate initial pressure head on each element
1631  // Assemble vector
1632  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_REVERSE);
1633  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_REVERSE);
1634  CHKERR VecAssemblyBegin(D1);
1635  CHKERR VecAssemblyEnd(D1);
1636  // Calculate norm
1637  double norm2D1;
1638  CHKERR VecNorm(D1, NORM_2, &norm2D1);
1639  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1640  PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
1642  }
1643 
1644  /**
1645  * \brief solve problem
1646  * @return error code
1647  */
1648  MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
1650  if (set_initial_pc) {
1651  // Set initial head
1652  CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
1653  }
1654 
1655  // Initiate vector from data on the mesh
1656  CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
1657 
1658  // Create time solver
1659  TS ts;
1660  CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
1661  // Use backward Euler method
1662  CHKERR TSSetType(ts, TSBEULER);
1663  // Set final time
1664  double ftime = 1;
1665  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
1666  // Setup solver from commabd line
1667  CHKERR TSSetFromOptions(ts);
1668  // Set DM to TS
1669  CHKERR TSSetDM(ts, dM);
1670 #if PETSC_VERSION_GE(3, 7, 0)
1671  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1672 #endif
1673  // Set-up monitor
1674  TsCtx *ts_ctx;
1675  DMMoFEMGetTsCtx(dM, &ts_ctx);
1676  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1677 
1678  // This add SNES monitor, to show error by fields. It is dirty trick
1679  // to add monitor, so code is hidden from doxygen
1680  CHKERR TSSetSolution(ts, D);
1681  CHKERR TSSetUp(ts);
1682  SNES snes;
1683  CHKERR TSGetSNES(ts, &snes);
1684 
1685 #if PETSC_VERSION_GE(3, 7, 0)
1686  {
1687  PetscViewerAndFormat *vf;
1688  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1689  PETSC_VIEWER_DEFAULT, &vf);
1690  CHKERR SNESMonitorSet(
1691  snes,
1692  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1693  void *))SNESMonitorFields,
1694  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1695  }
1696 #else
1697  {
1698  CHKERR SNESMonitorSet(snes,
1699  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1700  void *))SNESMonitorFields,
1701  0, 0);
1702  }
1703 #endif
1704 
1705  CHKERR TSSolve(ts, D);
1706 
1707  // Get statistic form TS and print it
1708  CHKERR TSGetTime(ts, &ftime);
1709  PetscInt steps, snesfails, rejects, nonlinits, linits;
1710  CHKERR TSGetTimeStepNumber(ts, &steps);
1711  CHKERR TSGetSNESFailures(ts, &snesfails);
1712  CHKERR TSGetStepRejections(ts, &rejects);
1713  CHKERR TSGetSNESIterations(ts, &nonlinits);
1714  CHKERR TSGetKSPIterations(ts, &linits);
1715  PetscPrintf(PETSC_COMM_WORLD,
1716  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
1717  "%D, linits %D\n",
1718  steps, rejects, snesfails, ftime, nonlinits, linits);
1719 
1721  }
1722 };
1723 
1724 } // namespace MixTransport
1725 
1726 #endif // __UNSATURATD_FLOW_HPP__
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Integrate boundary condition.
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
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
MoFEMErrorCode postProcess()
function is run at the end of loop
Pre-peprocessing Set head pressute rate and get inital essential boundary conditions.
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
Vec D1
Vector with inital head capillary pressure.
double diffK
Derivative of hydraulic conductivity [L/s * L^2/F].
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
boost::shared_ptr< ForcesAndSourcesCore > feVolInitialPc
Calculate inital boundary conditions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
add fields
virtual void printMatParameters(const int id, const std::string &prefix) const =0
double diffC
Derivative of capacity [S^2/L^2 * L^2/F ].
virtual MoFEMErrorCode calTheta()
Force scale operator for reading two columns.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
field with continuous normal traction
Definition: definitions.h:173
virtual moab::Interface & get_moab()=0
Vec ghostFlux
Ghost Vector of integrated fluxes.
MoFEMErrorCode solveProblem(bool set_initial_pc=true)
solve problem
MoFEMErrorCode calculateInitialPc()
Calculate inital pressure head distribution.
double h_t
rate of hydraulic head
double C
Capacity [S^2/L^2].
VectorDouble valuesAtGaussPts
values at integration points on element
static double ePsilon0
Regularization parameter.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:434
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:414
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
Implementation of operators, problem and finite elements for unsaturated flow.
Evaluate boundary condition at the boundary.
OpIntegrateFluxes(UnsaturatedFlowElement &ctx, const std::string fluxes_name)
Constructor.
boost::shared_ptr< ForcesAndSourcesCore > feFaceRhs
Face element apply natural bc.
MoFEMErrorCode createMatrices()
Create vectors and matrices.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
virtual MoFEMErrorCode getMaterial(const EntityHandle ent, int &block_id) const
For given element handle get material block Id.
std::map< int, boost::shared_ptr< GenericMaterial > > MaterialsDoubleMap
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
OpDivTauU_HdivL2(UnsaturatedFlowElement &ctx, const std::string &flux_name_col, const std::string &val_name_row)
Constructor.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:190
boost::shared_ptr< ForcesAndSourcesCore > feFaceBc
Elemnet to calculate essential bc.
Range tEts
Elements with this material.
boost::shared_ptr< ForcesAndSourcesCore > fluxIntegrate
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.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
boost::shared_ptr< FEMethod > tsMonitor
Set integration rule to boundary elements.
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:761
MaterialsDoubleMap dMatMap
materials database
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
virtual MoFEMErrorCode calDiffC()
VectorDouble divergenceAtGaussPts
divergence at integration points on element
boost::shared_ptr< VectorDouble > headRateAtGaussPts
Vector keeps head rate.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
boost::shared_ptr< ForcesAndSourcesCore > feVolRhs
Assemble residual vector.
virtual MoFEMErrorCode calDiffK()
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Class storing information about boundary condition.
Set integration rule to volume elements.
Mix transport problemNote to solve this system you need to use direct solver or proper preconditioner...
PetscReal ts_t
time
MatrixDouble fluxesAtGaussPts
fluxes at integration points on element
virtual double initalPcEval() const =0
Initialize head.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
preProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
MoFEMErrorCode setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule=VolRule(), ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule())
Create finite element instances.
MoFEMErrorCode destroyMatrices()
Delete matrices and vector when no longer needed.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Integrate boundary condition.
const int fRequency
bool sYmm
If true assume that matrix is symmetric structure.
virtual int get_comm_rank() const =0
postProcessVol(UnsaturatedFlowElement &ctx, boost::shared_ptr< ForcesAndSourcesCore > &fe_ptr)
moab::Interface & postProcMesh
BcMap bcValueMap
Store boundary condition for head capillary pressure.
PetscReal ts_a
shift for U_tt (see PETSc Time Solver)
Post proces method for volume element Assemble vectors and matrices and apply essential boundary cond...
OpEvaluateInitiallHead(UnsaturatedFlowElement &ctx, const std::string &val_name)
Vec ts_F
residual vector
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble matrix.
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
TS ts
time solver
static double ePsilon1
Regularization parameter.
DataForcesAndSourcesCore::EntData EntData
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.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
virtual MoFEMErrorCode calSe()
auto getFTensor0IntegrationWeight()
Get integration weights.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
double K
Hydraulic conductivity [L/s].
Post processing.
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
map< int, boost::shared_ptr< BcData > > BcMap
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:109
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Do calculations.
double Se
Effective saturation.
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
boost::shared_ptr< MethodForForceScaling > scaleMethodFlux
Method scaling fluxes.
boost::function< double(const double x, const double y, const double z)> hookFun
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
OpRhsBcOnValues(UnsaturatedFlowElement &ctx, const std::string fluxes_name, boost::shared_ptr< MethodForForceScaling > &value_scale)
Constructor.
double sCale
Scale time dependent eq.
boost::shared_ptr< MethodForForceScaling > scaleMethodValue
Method scaling values.
MoFEMErrorCode buildProblem(BitRefLevel ref_level=BitRefLevel().set(0))
Build problem.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode calC()
UnsaturatedFlowElement & cTx
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const int gg, const double x, const double y, const double z, double &value)
Get value on boundary.
virtual MoFEMErrorCode calK()
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:990
boost::shared_ptr< ForcesAndSourcesCore > fePtr
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
OpVDivSigma_L2Hdiv(UnsaturatedFlowElement &ctx, const std::string &val_name_row, const std::string &flux_name_col)
Constructor.
OpPostProcMaterial(UnsaturatedFlowElement &ctx, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name)
boost::shared_ptr< ForcesAndSourcesCore > feVolLhs
Assemble tangent matrix.
DM dM
Discrete manager for unsaturated flow problem.
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
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:708
Mat ts_B
Preconditioner for ts_A.
boost::shared_ptr< ForcesAndSourcesCore > fePtr
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
UnsaturatedFlowElement(MoFEM::Interface &m_field)
virtual MPI_Comm & get_comm() const =0
int operator()(int, int, int p_data) const
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
boost::shared_ptr< MethodForForceScaling > valueScale
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble matrix.
Generic material model for unsaturated water transport.
MonitorPostProc(UnsaturatedFlowElement &ctx, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc, boost::shared_ptr< ForcesAndSourcesCore > flux_Integrate, const int frequency)
UnsaturatedFlowElement & cTx
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:971
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Do calculations.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
MoFEMErrorCode getDivergenceOfHDivBaseFunctions(int side, EntityType type, DataForcesAndSourcesCore::EntData &data, int gg, VectorDouble &div)
Get divergence of base functions at integration point.
OpTauDotSigma_HdivHdiv(UnsaturatedFlowElement &ctx, const std::string flux_name)
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
int operator()(int p_row, int p_col, int p_data) const
std::vector< EntityHandle > & mapGaussPts
OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
MoFEMErrorCode calculateEssentialBc()
Calculate boundary conditions for fluxes.
field with C-1 continuity
Definition: definitions.h:174