v0.9.1
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  // Get material parameters
359  int nb_gauss_pts = data.getN().size1();
360  for (int gg = 0; gg != nb_gauss_pts; gg++) {
361  // Get divergence
363  const double alpha = t_w * vol;
364  block_data->h = t_h;
365  block_data->x = t_coords(0);
366  block_data->y = t_coords(1);
367  block_data->z = t_coords(2);
368  CHKERR block_data->calK();
369  const double K = block_data->K;
370  const double scaleZ = block_data->scaleZ;
371  const double z = t_coords(2); /// z-coordinate at Gauss pt
372  // Calculate pressure gradient
373  noalias(nF) -= alpha * (t_h - z * scaleZ) * divVec;
374  // Calculate presure gradient from flux
375  FTensor::Tensor0<double *> t_nf(&*nF.begin());
376  for (int rr = 0; rr != nb_dofs; rr++) {
377  t_nf += alpha * (1 / K) * (t_n_hdiv(i) * t_flux(i));
378  ++t_n_hdiv; // move to next base function
379  ++t_nf; // move to next element in vector
380  }
381  ++t_h; // move to next integration point
382  ++t_flux;
383  ++t_coords;
384  ++t_w;
385  }
386  // Assemble residual
387  CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs,
388  &*data.getIndices().begin(), &*nF.begin(),
389  ADD_VALUES);
391  }
392  };
393 
394  /**
395  * Assemble mass residual
396  */
399 
401 
402  OpResidualMass(UnsaturatedFlowElement &ctx, const std::string &val_name)
404  val_name, UserDataOperator::OPROW),
405  cTx(ctx) {}
406 
408 
409  MoFEMErrorCode doWork(int side, EntityType type,
412  const int nb_dofs = data.getIndices().size();
413  if (nb_dofs == 0)
415  // Resize local element vector
416  nF.resize(nb_dofs, false);
417  nF.clear();
418  // Get EntityHandle of the finite element
419  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
420  // Get material block id
421  int block_id;
422  CHKERR cTx.getMaterial(fe_ent, block_id);
423  // Get material block
424  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
425  // Get pressure
427  // Get pressure rate
429  // Flux divergence
430  auto t_div_flux = getFTensor0FromVec(cTx.divergenceAtGaussPts);
431  // Get integration weight
432  auto t_w = getFTensor0IntegrationWeight();
433  // Coords at integration points
434  auto t_coords = getFTensor1CoordsAtGaussPts();
435  // Scale eq.
436  const double scale = block_data->sCale;
437  // Get volume
438  const double vol = getVolume();
439  // Get number of integration points
440  int nb_gauss_pts = data.getN().size1();
441  for (int gg = 0; gg != nb_gauss_pts; gg++) {
442  const double alpha = t_w * vol * scale;
443  block_data->h = t_h;
444  block_data->x = t_coords(0);
445  block_data->y = t_coords(1);
446  block_data->z = t_coords(2);
447  CHKERR block_data->calC();
448  const double C = block_data->C;
449  // Calculate flux conservation
450  noalias(nF) += (alpha * (t_div_flux + C * t_h_t)) * data.getN(gg);
451  ++t_h;
452  ++t_h_t;
453  ++t_div_flux;
454  ++t_coords;
455  ++t_w;
456  }
457  // Assemble local vector
458  Vec f = getFEMethod()->ts_F;
459  CHKERR VecSetValues(f, nb_dofs, &*data.getIndices().begin(), &*nF.begin(),
460  ADD_VALUES);
462  }
463  };
464 
467 
469 
471  const std::string flux_name)
473  flux_name, flux_name, UserDataOperator::OPROWCOL),
474  cTx(ctx) {
475  sYmm = true;
476  }
477 
479 
481 
482  /**
483  * \brief Assemble matrix
484  * @param row_side local index of row entity on element
485  * @param col_side local index of col entity on element
486  * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
487  * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
488  * @param row_data data for row
489  * @param col_data data for col
490  * @return error code
491  */
492  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
493  EntityType col_type,
497  const int nb_row = row_data.getIndices().size();
498  const int nb_col = col_data.getIndices().size();
499  if (nb_row == 0)
501  if (nb_col == 0)
503  nN.resize(nb_row, nb_col, false);
504  nN.clear();
505  // Get EntityHandle of the finite element
506  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
507  // Get material block id
508  int block_id;
509  CHKERR cTx.getMaterial(fe_ent, block_id);
510  // Get material block
511  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
512  // Get pressure
514  // Coords at integration points
515  auto t_coords = getFTensor1CoordsAtGaussPts();
516  // Get base functions
517  auto t_n_hdiv_row = row_data.getFTensor1N<3>();
518  // Get integration weight
519  auto t_w = getFTensor0IntegrationWeight();
520  // Get volume
521  const double vol = getVolume();
522  int nb_gauss_pts = row_data.getN().size1();
523  for (int gg = 0; gg != nb_gauss_pts; gg++) {
524  block_data->h = t_h;
525  block_data->x = t_coords(0);
526  block_data->y = t_coords(1);
527  block_data->z = t_coords(2);
528  CHKERR block_data->calK();
529  const double K = block_data->K;
530  // get integration weight and multiply by element volume
531  const double alpha = t_w * vol;
532  const double beta = alpha * (1 / K);
533  FTensor::Tensor0<double *> t_a(&*nN.data().begin());
534  for (int kk = 0; kk != nb_row; kk++) {
535  auto t_n_hdiv_col = col_data.getFTensor1N<3>(gg, 0);
536  for (int ll = 0; ll != nb_col; ll++) {
537  t_a += beta * (t_n_hdiv_row(j) * t_n_hdiv_col(j));
538  ++t_n_hdiv_col;
539  ++t_a;
540  }
541  ++t_n_hdiv_row;
542  }
543  ++t_coords;
544  ++t_h;
545  ++t_w;
546  }
547  Mat a = getFEMethod()->ts_B;
548  CHKERR MatSetValues(a, nb_row, &*row_data.getIndices().begin(), nb_col,
549  &*col_data.getIndices().begin(), &*nN.data().begin(),
550  ADD_VALUES);
551  // matrix is symmetric, assemble other part
552  if (row_side != col_side || row_type != col_type) {
553  transNN.resize(nb_col, nb_row);
554  noalias(transNN) = trans(nN);
555  CHKERR MatSetValues(a, nb_col, &*col_data.getIndices().begin(), nb_row,
556  &*row_data.getIndices().begin(),
557  &*transNN.data().begin(), ADD_VALUES);
558  }
560  }
561  };
562 
563  struct OpVU_L2L2
565 
567 
568  OpVU_L2L2(UnsaturatedFlowElement &ctx, const std::string value_name)
570  value_name, value_name, UserDataOperator::OPROWCOL),
571  cTx(ctx) {
572  sYmm = true;
573  }
574 
576 
577  /**
578  * \brief Assemble matrix
579  * @param row_side local index of row entity on element
580  * @param col_side local index of col entity on element
581  * @param row_type type of row entity, f.e. MBVERTEX, MBEDGE, or MBTET
582  * @param col_type type of col entity, f.e. MBVERTEX, MBEDGE, or MBTET
583  * @param row_data data for row
584  * @param col_data data for col
585  * @return error code
586  */
587  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
588  EntityType col_type,
592  int nb_row = row_data.getIndices().size();
593  int nb_col = col_data.getIndices().size();
594  if (nb_row == 0)
596  if (nb_col == 0)
598  nN.resize(nb_row, nb_col, false);
599  nN.clear();
600  // Get EntityHandle of the finite element
601  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
602  // Get material block id
603  int block_id;
604  CHKERR cTx.getMaterial(fe_ent, block_id);
605  // Get material block
606  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
607  // Get pressure
609  // // Get pressure
610  // auto t_flux_residual = getFTensor0FromVec(*cTx.resAtGaussPts);
611  // Get pressure rate
613  // Get integration weight
614  auto t_w = getFTensor0IntegrationWeight();
615  // Coords at integration points
616  auto t_coords = getFTensor1CoordsAtGaussPts();
617  // Scale eq.
618  const double scale = block_data->sCale;
619  // Time step factor
620  double ts_a = getFEMethod()->ts_a;
621  // get volume
622  const double vol = getVolume();
623  int nb_gauss_pts = row_data.getN().size1();
624  // get base functions on rows
625  auto t_n_row = row_data.getFTensor0N();
626  for (int gg = 0; gg != nb_gauss_pts; gg++) {
627  // get integration weight and multiply by element volume
628  double alpha = t_w * vol * scale;
629  // evaluate material model at integration points
630  // to calculate capacity and tangent of capacity term
631  block_data->h = t_h;
632  block_data->h_t = t_h_t;
633  block_data->x = t_coords(0);
634  block_data->y = t_coords(1);
635  block_data->z = t_coords(2);
636  CHKERR block_data->calC();
637  CHKERR block_data->calDiffC();
638  const double C = block_data->C;
639  const double diffC = block_data->diffC;
640  // assemble local entity tangent matrix
642  &*nN.data().begin());
643  // iterate base functions on rows
644  for (int kk = 0; kk != nb_row; kk++) {
645  // get first base function on column at integration point gg
646  auto t_n_col = col_data.getFTensor0N(gg, 0);
647  // iterate base functions on columns
648  for (int ll = 0; ll != nb_col; ll++) {
649  // assemble elements of local matrix
650  t_a += (alpha * (C * ts_a + diffC * t_h_t)) * t_n_row * t_n_col;
651  ++t_n_col; // move to next base function on column
652  ++t_a; // move to next element in local tangent matrix
653  }
654  ++t_n_row; // move to next base function on row
655  }
656  ++t_w; // move to next integration weight
657  ++t_coords; // move to next coordinate at integration point
658  ++t_h; // move to next capillary head at integration point
659  // ++t_flux_residual;
660  ++t_h_t; // move to next capillary head rate at integration point
661  }
662  Mat a = getFEMethod()->ts_B;
663  CHKERR MatSetValues(a, nb_row, &row_data.getIndices()[0], nb_col,
664  &col_data.getIndices()[0], &*nN.data().begin(),
665  ADD_VALUES);
667  }
668  };
669 
672 
674 
675  /**
676  * \brief Constructor
677  */
679  const std::string &val_name_row,
680  const std::string &flux_name_col)
682  val_name_row, flux_name_col, UserDataOperator::OPROWCOL, false),
683  cTx(ctx) {}
684 
687 
688  /**
689  * \brief Do calculations
690  * @param row_side local index of entity on row
691  * @param col_side local index of entity on column
692  * @param row_type type of row entity
693  * @param col_type type of col entity
694  * @param row_data row data structure carrying information about base
695  * functions, DOFs indices, etc.
696  * @param col_data column data structure carrying information about base
697  * functions, DOFs indices, etc.
698  * @return error code
699  */
700  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
701  EntityType col_type,
705  int nb_row = row_data.getFieldData().size();
706  int nb_col = col_data.getFieldData().size();
707  if (nb_row == 0)
709  if (nb_col == 0)
711  // Get EntityHandle of the finite element
712  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
713  // Get material block id
714  int block_id;
715  CHKERR cTx.getMaterial(fe_ent, block_id);
716  // Get material block
717  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
718  nN.resize(nb_row, nb_col, false);
719  divVec.resize(nb_col, false);
720  nN.clear();
721  // Scale eq.
722  const double scale = block_data->sCale;
723  int nb_gauss_pts = row_data.getN().size1();
724  for (int gg = 0; gg < nb_gauss_pts; gg++) {
725  double alpha = getGaussPts()(3, gg) * getVolume() * scale;
726  CHKERR getDivergenceOfHDivBaseFunctions(col_side, col_type, col_data,
727  gg, divVec);
728  noalias(nN) += alpha * outer_prod(row_data.getN(gg), divVec);
729  }
730  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
731  &row_data.getIndices()[0], nb_col,
732  &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
734  }
735  };
736 
739 
741 
742  /**
743  * \brief Constructor
744  */
746  const std::string &flux_name_col,
747  const std::string &val_name_row)
749  flux_name_col, val_name_row, UserDataOperator::OPROWCOL, false),
750  cTx(ctx) {}
751 
755 
756  /**
757  * \brief Do calculations
758  * @param row_side local index of entity on row
759  * @param col_side local index of entity on column
760  * @param row_type type of row entity
761  * @param col_type type of col entity
762  * @param row_data row data structure carrying information about base
763  * functions, DOFs indices, etc.
764  * @param col_data column data structure carrying information about base
765  * functions, DOFs indices, etc.
766  * @return error code
767  */
768  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
769  EntityType col_type,
773  int nb_row = row_data.getFieldData().size();
774  int nb_col = col_data.getFieldData().size();
775  if (nb_row == 0)
777  if (nb_col == 0)
779  nN.resize(nb_row, nb_col, false);
780  divVec.resize(nb_row, false);
781  nN.clear();
782  // Get EntityHandle of the finite element
783  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
784  // Get material block id
785  int block_id;
786  CHKERR cTx.getMaterial(fe_ent, block_id);
787  // Get material block
788  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
789  // Get pressure
791  // Get flux
792  auto t_flux = getFTensor1FromMat<3>(cTx.fluxesAtGaussPts);
793  // Coords at integration points
794  auto t_coords = getFTensor1CoordsAtGaussPts();
795  // Get integration weight
796  auto t_w = getFTensor0IntegrationWeight();
797  // Get base function
798  auto t_n_hdiv_row = row_data.getFTensor1N<3>();
799  // Get volume
800  double vol = getVolume();
801  int nb_gauss_pts = row_data.getN().size1();
802  for (int gg = 0; gg != nb_gauss_pts; gg++) {
803  block_data->h = t_h;
804  block_data->x = t_coords(0);
805  block_data->y = t_coords(1);
806  block_data->z = t_coords(2);
807  CHKERR block_data->calK();
808  CHKERR block_data->calDiffK();
809  const double K = block_data->K;
810  // const double z = t_coords(2);
811  const double KK = K * K;
812  const double diffK = block_data->diffK;
813  double alpha = t_w * vol;
814  CHKERR getDivergenceOfHDivBaseFunctions(row_side, row_type, row_data,
815  gg, divVec);
816  noalias(nN) -= alpha * outer_prod(divVec, col_data.getN(gg));
817  FTensor::Tensor0<double *> t_a(&*nN.data().begin());
818  for (int rr = 0; rr != nb_row; rr++) {
819  double beta = alpha * (-diffK / KK) * (t_n_hdiv_row(i) * t_flux(i));
820  auto t_n_col = col_data.getFTensor0N(gg, 0);
821  for (int cc = 0; cc != nb_col; cc++) {
822  t_a += beta * t_n_col;
823  ++t_n_col;
824  ++t_a;
825  }
826  ++t_n_hdiv_row;
827  }
828  ++t_w;
829  ++t_coords;
830  ++t_h;
831  ++t_flux;
832  }
833  CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
834  &row_data.getIndices()[0], nb_col,
835  &col_data.getIndices()[0], &nN(0, 0), ADD_VALUES);
837  }
838  };
839 
844  const std::string &val_name)
846  val_name, UserDataOperator::OPROW),
847  cTx(ctx) {}
848 
851 
852  MoFEMErrorCode doWork(int side, EntityType type,
855  if (data.getFieldData().size() == 0)
857  int nb_dofs = data.getFieldData().size();
858  int nb_gauss_pts = data.getN().size1();
859  if (nb_dofs != static_cast<int>(data.getN().size2())) {
860  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
861  "wrong number of dofs");
862  }
863  nN.resize(nb_dofs, nb_dofs, false);
864  nF.resize(nb_dofs, false);
865  nN.clear();
866  nF.clear();
867 
868  // Get EntityHandle of the finite element
869  EntityHandle fe_ent = getFEEntityHandle();
870  // Get material block id
871  int block_id;
872  CHKERR cTx.getMaterial(fe_ent, block_id);
873  // Get material block
874  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
875 
876  // loop over integration points
877  for (int gg = 0; gg < nb_gauss_pts; gg++) {
878  // get coordinates at integration point
879  block_data->x = getCoordsAtGaussPts()(gg, 0);
880  block_data->y = getCoordsAtGaussPts()(gg, 1);
881  block_data->z = getCoordsAtGaussPts()(gg, 2);
882  // get weight for integration rule
883  double alpha = getGaussPts()(2, gg) * getVolume();
884  nN += alpha * outer_prod(data.getN(gg), data.getN(gg));
885  nF += alpha * block_data->initialPcEval() * data.getN(gg);
886  }
887 
888  // factor matrix
890  // solve local problem
891  cholesky_solve(nN, nF, ublas::lower());
892 
893  // set solution to vector
894  CHKERR VecSetValues(cTx.D1, nb_dofs, &*data.getIndices().begin(),
895  &*nF.begin(), INSERT_VALUES);
896 
898  }
899  };
900 
904 
905  /**
906  * \brief Constructor
907  */
909  const std::string fluxes_name)
911  fluxes_name, UserDataOperator::OPROW),
912  cTx(ctx) {}
913 
915 
916  /**
917  * \brief Integrate boundary condition
918  * @param side local index of entity
919  * @param type type of entity
920  * @param data data on entity
921  * @return error code
922  */
923  MoFEMErrorCode doWork(int side, EntityType type,
926  int nb_dofs = data.getFieldData().size();
927  if (nb_dofs == 0)
929  // Get base function
930  auto t_n_hdiv = data.getFTensor1N<3>();
931  // get normal of face
932  auto t_normal = getFTensor1Normal();
933  // Integration weight
934  auto t_w = getFTensor0IntegrationWeight();
935  double flux_on_entity = 0;
936  int nb_gauss_pts = data.getN().size1();
937  for (int gg = 0; gg < nb_gauss_pts; gg++) {
938  auto t_data = data.getFTensor0FieldData();
939  for (int rr = 0; rr != nb_dofs; rr++) {
940  flux_on_entity -= (0.5 * t_data * t_w) * (t_n_hdiv(i) * t_normal(i));
941  ++t_n_hdiv;
942  ++t_data;
943  }
944  ++t_w;
945  }
946  CHKERR VecSetValue(cTx.ghostFlux, 0, flux_on_entity, ADD_VALUES);
948  }
949  };
950 
951  /**
952  * Operator used to post-process results for unsaturated infiltration problem.
953  * Operator should with element for post-processing results, i.e.
954  * PostProcVolumeOnRefinedMesh
955  */
958 
961  std::vector<EntityHandle> &mapGaussPts;
962 
964  moab::Interface &post_proc_mesh,
965  std::vector<EntityHandle> &map_gauss_pts,
966  const std::string field_name)
969  cTx(ctx), postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts) {}
970 
971  MoFEMErrorCode doWork(int side, EntityType type,
974  int nb_dofs = data.getFieldData().size();
975  if (nb_dofs == 0)
977 
978  // Get EntityHandle of the finite element
979  EntityHandle fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
980  // Get material block id
981  int block_id;
982  CHKERR cTx.getMaterial(fe_ent, block_id);
983  // Get material block
984  boost::shared_ptr<GenericMaterial> &block_data = cTx.dMatMap.at(block_id);
985 
986  // Set bloc Id
987  Tag th_id;
988  int def_block_id = -1;
989  CHKERR postProcMesh.tag_get_handle("BLOCK_ID", 1, MB_TYPE_INTEGER, th_id,
990  MB_TAG_CREAT | MB_TAG_SPARSE,
991  &def_block_id);
992  CHKERR postProcMesh.tag_set_data(th_id, &fe_ent, 1, &block_id);
993 
994  // Create mesh tag. Tags are created on post-processing mesh and
995  // visable in post-processor, e.g. Paraview
996  double zero = 0;
997  Tag th_theta;
998  CHKERR postProcMesh.tag_get_handle("THETA", 1, MB_TYPE_DOUBLE, th_theta,
999  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1000  Tag th_se;
1001  CHKERR postProcMesh.tag_get_handle("Se", 1, MB_TYPE_DOUBLE, th_se,
1002  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1003  // Tag th_ks;
1004  // CHKERR postProcMesh.tag_get_handle(
1005  // "Ks",1,MB_TYPE_DOUBLE,th_ks,
1006  // MB_TAG_CREAT|MB_TAG_SPARSE,&zero
1007  // );
1008  // CHKERR postProcMesh.tag_set_data(th_ks,&fe_ent,1,&block_data->Ks);
1009  Tag th_k;
1010  CHKERR postProcMesh.tag_get_handle("K", 1, MB_TYPE_DOUBLE, th_k,
1011  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1012  Tag th_c;
1013  CHKERR postProcMesh.tag_get_handle("C", 1, MB_TYPE_DOUBLE, th_c,
1014  MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1015 
1016  // Get pressure at integration points
1018  // Coords at integration points
1019  auto t_coords = getFTensor1CoordsAtGaussPts();
1020 
1021  int nb_gauss_pts = data.getN().size1();
1022  for (int gg = 0; gg < nb_gauss_pts; gg++) {
1023  block_data->h = t_h;
1024  block_data->x = t_coords(0);
1025  block_data->y = t_coords(1);
1026  block_data->z = t_coords(2);
1027  // Calculate theta (water content) and save it on mesh tags
1028  CHKERR block_data->calTheta();
1029  double theta = block_data->tHeta;
1030  CHKERR postProcMesh.tag_set_data(th_theta, &mapGaussPts[gg], 1, &theta);
1031  CHKERR block_data->calSe();
1032  // Calculate Se (effective saturation and save it on the mesh tags)
1033  double Se = block_data->Se;
1034  CHKERR postProcMesh.tag_set_data(th_se, &mapGaussPts[gg], 1, &Se);
1035  // Calculate K (hydraulic conductivity) and save it on the mesh tags
1036  CHKERR block_data->calK();
1037  double K = block_data->K;
1038  CHKERR postProcMesh.tag_set_data(th_k, &mapGaussPts[gg], 1, &K);
1039  // Calculate water capacity and save it on the mesh tags
1040  CHKERR block_data->calC();
1041  double C = block_data->C;
1042  CHKERR postProcMesh.tag_set_data(th_c, &mapGaussPts[gg], 1, &C);
1043  ++t_h;
1044  ++t_coords;
1045  }
1046 
1048  }
1049  };
1050 
1051  /**
1052  * Finite element implementation called by TS monitor. Element calls other
1053  * finite elements to evaluate material properties and save results on the
1054  * mesh.
1055  *
1056  * \note Element overloaded only FEMethod::postProcess methos where other
1057  * elements are called.
1058  */
1060 
1062  boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProc;
1063  boost::shared_ptr<ForcesAndSourcesCore> fluxIntegrate;
1064 
1065  const int fRequency;
1066 
1068  boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc,
1069  boost::shared_ptr<ForcesAndSourcesCore> flux_Integrate,
1070  const int frequency)
1071  : cTx(ctx), postProc(post_proc), fluxIntegrate(flux_Integrate),
1072  fRequency(frequency) {}
1073 
1077  }
1078 
1082  }
1083 
1086 
1087  // Get time step
1088  int step;
1089  CHKERR TSGetTimeStepNumber(ts, &step);
1090 
1091  if ((step) % fRequency == 0) {
1092  // Post-process results and save in the file
1093  PetscPrintf(PETSC_COMM_WORLD, "Output results %d - %d\n", step,
1094  fRequency);
1096  CHKERR postProc->writeFile(
1097  string("out_") + boost::lexical_cast<std::string>(step) + ".h5m");
1098  }
1099 
1100  // Integrate fluxes on faces where pressure head is applied
1101  CHKERR VecZeroEntries(cTx.ghostFlux);
1102  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1103  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1104  // Run finite element to integrate fluxes
1106  CHKERR VecAssemblyBegin(cTx.ghostFlux);
1107  CHKERR VecAssemblyEnd(cTx.ghostFlux);
1108  // accumulate errors from processors
1109  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1110  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, ADD_VALUES, SCATTER_REVERSE);
1111  // scatter errors to all processors
1112  CHKERR VecGhostUpdateBegin(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1113  CHKERR VecGhostUpdateEnd(cTx.ghostFlux, INSERT_VALUES, SCATTER_FORWARD);
1114  double *ghost_flux;
1115  CHKERR VecGetArray(cTx.ghostFlux, &ghost_flux);
1116  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Flux at time %6.4g %6.4g\n", ts_t,
1117  ghost_flux[0]);
1118  CHKERR VecRestoreArray(cTx.ghostFlux, &ghost_flux);
1119 
1121  }
1122  };
1123 
1124  /// \brief add fields
1125  MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes,
1126  const int order) {
1128  // Fields
1131  CHKERR mField.add_field(values + "_t", L2, AINSWORTH_LEGENDRE_BASE, 1);
1132  // CHKERR mField.add_field(fluxes+"_residual",L2,AINSWORTH_LEGENDRE_BASE,1);
1133 
1134  // meshset consisting all entities in mesh
1135  EntityHandle root_set = mField.get_moab().get_root_set();
1136  // add entities to field
1137 
1139  if (it->getName().compare(0, 4, "SOIL") != 0)
1140  continue;
1141  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1142  MBTET, fluxes);
1143  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1144  MBTET, values);
1145  CHKERR mField.add_ents_to_field_by_type(dMatMap[it->getMeshsetId()]->tEts,
1146  MBTET, values + "_t");
1147  // CHKERR mField.add_ents_to_field_by_type(
1148  // dMatMap[it->getMeshsetId()]->tEts,MBTET,fluxes+"_residual"
1149  // );
1150  }
1151 
1152  CHKERR mField.set_field_order(root_set, MBTET, fluxes, order + 1);
1153  CHKERR mField.set_field_order(root_set, MBTRI, fluxes, order + 1);
1154  CHKERR mField.set_field_order(root_set, MBTET, values, order);
1155  CHKERR mField.set_field_order(root_set, MBTET, values + "_t", order);
1156  // CHKERR mField.set_field_order(root_set,MBTET,fluxes+"_residual",order);
1158  }
1159 
1160  /// \brief add finite elements
1161  MoFEMErrorCode addFiniteElements(const std::string &fluxes_name,
1162  const std::string &values_name) {
1164 
1165  // Define element "MIX". Note that this element will work with fluxes_name
1166  // and values_name. This reflect bilinear form for the problem
1175  values_name + "_t");
1176  // CHKERR
1177  // mField.modify_finite_element_add_field_data("MIX",fluxes_name+"_residual");
1178 
1180  if (it->getName().compare(0, 4, "SOIL") != 0)
1181  continue;
1183  dMatMap[it->getMeshsetId()]->tEts, MBTET, "MIX");
1184  }
1185 
1186  // Define element to integrate natural boundary conditions, i.e. set values.
1187  CHKERR mField.add_finite_element("MIX_BCVALUE", MF_ZERO);
1189  fluxes_name);
1191  fluxes_name);
1193  fluxes_name);
1195  values_name);
1196 
1198  if (it->getName().compare(0, 4, "HEAD") != 0)
1199  continue;
1201  bcValueMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCVALUE");
1202  }
1203 
1204  // Define element to apply essential boundary conditions.
1205  CHKERR mField.add_finite_element("MIX_BCFLUX", MF_ZERO);
1207  fluxes_name);
1209  fluxes_name);
1211  fluxes_name);
1213  values_name);
1214 
1216  if (it->getName().compare(0, 4, "FLUX") != 0)
1217  continue;
1219  bcFluxMap[it->getMeshsetId()]->eNts, MBTRI, "MIX_BCFLUX");
1220  }
1221 
1223  }
1224 
1225  /**
1226  * \brief Build problem
1227  * @param ref_level mesh refinement on which mesh problem you like to built.
1228  * @return error code
1229  */
1230  MoFEMErrorCode buildProblem(Range zero_flux_ents,
1231  BitRefLevel ref_level = BitRefLevel().set(0)) {
1233 
1234  // Build fields
1236  // Build finite elements
1238  CHKERR mField.build_finite_elements("MIX_BCFLUX");
1239  CHKERR mField.build_finite_elements("MIX_BCVALUE");
1240  // Build adjacencies of degrees of freedom and elements
1241  CHKERR mField.build_adjacencies(ref_level);
1242 
1243  // create DM instance
1244  CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
1245  // setting that DM is type of DMMOFEM, i.e. MOFEM implementation manages DM
1246  CHKERR DMSetType(dM, "DMMOFEM");
1247  // mesh is portioned, each process keeps only part of problem
1248  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1249  // creates problem in DM
1250  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "MIX", ref_level);
1251  // discretised problem creates square matrix (that makes some optimizations)
1252  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1253  // set DM options from command line
1254  CHKERR DMSetFromOptions(dM);
1255  // add finite elements
1256  CHKERR DMMoFEMAddElement(dM, "MIX");
1257  CHKERR DMMoFEMAddElement(dM, "MIX_BCFLUX");
1258  CHKERR DMMoFEMAddElement(dM, "MIX_BCVALUE");
1259  // constructor data structures
1260  CHKERR DMSetUp(dM);
1261 
1262  // remove zero flux dofs
1263  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
1264  "MIX", "FLUXES", zero_flux_ents);
1265 
1266  PetscSection section;
1267  CHKERR mField.getInterface<ISManager>()->sectionCreate("MIX", &section);
1268  CHKERR DMSetDefaultSection(dM, section);
1269  CHKERR DMSetDefaultGlobalSection(dM, section);
1270  // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
1271  CHKERR PetscSectionDestroy(&section);
1272 
1274  }
1275 
1276  boost::shared_ptr<ForcesAndSourcesCore>
1277  feFaceBc; ///< Elemnet to calculate essential bc
1278  boost::shared_ptr<ForcesAndSourcesCore>
1279  feFaceRhs; ///< Face element apply natural bc
1280  boost::shared_ptr<ForcesAndSourcesCore>
1281  feVolInitialPc; ///< Calculate inital boundary conditions
1282  boost::shared_ptr<ForcesAndSourcesCore>
1283  feVolRhs; ///< Assemble residual vector
1284  boost::shared_ptr<ForcesAndSourcesCore> feVolLhs; ///< Assemble tangent matrix
1285  boost::shared_ptr<MethodForForceScaling>
1286  scaleMethodFlux; ///< Method scaling fluxes
1287  boost::shared_ptr<MethodForForceScaling>
1288  scaleMethodValue; ///< Method scaling values
1289  boost::shared_ptr<FEMethod> tsMonitor; ///< Element used by TS monitor to
1290  ///< postprocess results at time step
1291 
1292  boost::shared_ptr<VectorDouble>
1293  headRateAtGaussPts; ///< Vector keeps head rate
1294  // boost::shared_ptr<VectorDouble> resAtGaussPts; ///< Residual field
1295 
1296  /**
1297  * \brief Set integration rule to volume elements
1298  *
1299  */
1300  struct VolRule {
1301  int operator()(int, int, int p_data) const { return 2 * p_data + p_data; }
1302  };
1303  /**
1304  * \brief Set integration rule to boundary elements
1305  *
1306  */
1307  struct FaceRule {
1308  int operator()(int p_row, int p_col, int p_data) const {
1309  return 2 * p_data;
1310  }
1311  };
1312 
1313  std::vector<int> bcVecIds;
1315 
1316  /**
1317  * \brief Pre-peprocessing
1318  * Set head pressute rate and get inital essential boundary conditions
1319  */
1320  struct preProcessVol {
1322  boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1323  // std::string mArk;
1324 
1327  boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr /*,std::string mark*/
1328  )
1329  : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1332  // Update pressure rates
1333  CHKERR fePtr->mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1334  fePtr->problemPtr, "VALUES", string("VALUES") + "_t", ROW,
1335  fePtr->ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1336  switch (fePtr->ts_ctx) {
1337  case TSMethod::CTX_TSSETIFUNCTION:
1338  CHKERR VecSetOption(fePtr->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
1339  PETSC_TRUE);
1340  if (!cTx.bcIndices.empty()) {
1341  double scale;
1342  CHKERR cTx.scaleMethodFlux->getForceScale(fePtr->ts_t, scale);
1343  if (cTx.bcVecIds.size() != cTx.bcIndices.size()) {
1344  cTx.bcVecIds.insert(cTx.bcVecIds.begin(), cTx.bcIndices.begin(),
1345  cTx.bcIndices.end());
1346  cTx.bcVecVals.resize(cTx.bcVecIds.size(), false);
1347  cTx.vecValsOnBc.resize(cTx.bcVecIds.size(), false);
1348  }
1349  CHKERR VecGetValues(cTx.D0, cTx.bcVecIds.size(),
1350  &*cTx.bcVecIds.begin(), &*cTx.bcVecVals.begin());
1351  CHKERR VecGetValues(fePtr->ts_u, cTx.bcVecIds.size(),
1352  &*cTx.bcVecIds.begin(),
1353  &*cTx.vecValsOnBc.begin());
1354  cTx.bcVecVals *= scale;
1355  VectorDouble::iterator vit = cTx.bcVecVals.begin();
1356  const NumeredDofEntity *dof_ptr;
1357  for (std::vector<int>::iterator it = cTx.bcVecIds.begin();
1358  it != cTx.bcVecIds.end(); it++, vit++) {
1359  if (auto dof_ptr =
1360  fePtr->problemPtr->getColDofsByPetscGlobalDofIdx(*it)
1361  .lock())
1362  dof_ptr->getFieldData() = *vit;
1363  }
1364  } else {
1365  cTx.bcVecIds.resize(0);
1366  cTx.bcVecVals.resize(0);
1367  cTx.vecValsOnBc.resize(0);
1368  }
1369  break;
1370  default:
1371  // don nothing
1372  break;
1373  }
1375  }
1376  };
1377 
1378  /**
1379  * \brief Post proces method for volume element
1380  * Assemble vectors and matrices and apply essential boundary conditions
1381  */
1384  boost::shared_ptr<ForcesAndSourcesCore> fePtr;
1385  // std::string mArk;
1388  boost::shared_ptr<ForcesAndSourcesCore> &fe_ptr //,std::string mark
1389  )
1390  : cTx(ctx), fePtr(fe_ptr) /*,mArk(mark)*/ {}
1393  switch (fePtr->ts_ctx) {
1394  case TSMethod::CTX_TSSETIJACOBIAN: {
1395  CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1396  CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1397  // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1398  // std::string wait;
1399  // std::cin >> wait;
1400  CHKERR MatZeroRowsColumns(fePtr->ts_B, cTx.bcVecIds.size(),
1401  &*cTx.bcVecIds.begin(), 1, PETSC_NULL,
1402  PETSC_NULL);
1403  CHKERR MatAssemblyBegin(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1404  CHKERR MatAssemblyEnd(fePtr->ts_B, MAT_FINAL_ASSEMBLY);
1405  // MatView(fePtr->ts_B,PETSC_VIEWER_DRAW_WORLD);
1406  // std::string wait;
1407  // std::cin >> wait;
1408  } break;
1409  case TSMethod::CTX_TSSETIFUNCTION: {
1410  CHKERR VecAssemblyBegin(fePtr->ts_F);
1411  CHKERR VecAssemblyEnd(fePtr->ts_F);
1412  if (!cTx.bcVecIds.empty()) {
1414  CHKERR VecSetValues(fePtr->ts_F, cTx.bcVecIds.size(),
1415  &*cTx.bcVecIds.begin(), &*cTx.vecValsOnBc.begin(),
1416  INSERT_VALUES);
1417  }
1418  CHKERR VecAssemblyBegin(fePtr->ts_F);
1419  CHKERR VecAssemblyEnd(fePtr->ts_F);
1420  } break;
1421  default:
1422  // don nothing
1423  break;
1424  }
1426  }
1427  };
1428 
1429  /**
1430  * \brief Create finite element instances
1431  * @param vol_rule integration rule for volume element
1432  * @param face_rule integration rule for boundary element
1433  * @return error code
1434  */
1436  setFiniteElements(ForcesAndSourcesCore::RuleHookFun vol_rule = VolRule(),
1437  ForcesAndSourcesCore::RuleHookFun face_rule = FaceRule()) {
1439 
1440  // create finite element instances
1441  feFaceBc = boost::shared_ptr<ForcesAndSourcesCore>(
1443  feFaceRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1445  feVolInitialPc = boost::shared_ptr<ForcesAndSourcesCore>(
1447  feVolLhs = boost::shared_ptr<ForcesAndSourcesCore>(
1449  feVolRhs = boost::shared_ptr<ForcesAndSourcesCore>(
1451  // set integration rule to elements
1452  feFaceBc->getRuleHook = face_rule;
1453  feFaceRhs->getRuleHook = face_rule;
1454  feVolInitialPc->getRuleHook = vol_rule;
1455  feVolLhs->getRuleHook = vol_rule;
1456  feVolRhs->getRuleHook = vol_rule;
1457  // set function hook for finite element postprocessing stage
1458  feVolRhs->preProcessHook = preProcessVol(*this, feVolRhs);
1459  feVolLhs->preProcessHook = preProcessVol(*this, feVolLhs);
1460  feVolRhs->postProcessHook = postProcessVol(*this, feVolRhs);
1461  feVolLhs->postProcessHook = postProcessVol(*this, feVolLhs);
1462 
1463  // create method for setting history for fluxes on boundary
1464  scaleMethodFlux = boost::shared_ptr<MethodForForceScaling>(
1465  new TimeForceScale("-flux_history", false));
1466 
1467  // create method for setting history for presure heads on boundary
1468  scaleMethodValue = boost::shared_ptr<MethodForForceScaling>(
1469  new TimeForceScale("-head_history", false));
1470 
1471  // Set operator to calculate essential boundary conditions
1472  feFaceBc->getOpPtrVector().push_back(
1473  new OpEvaluateBcOnFluxes(*this, "FLUXES"));
1474 
1475  // Set operator to calculate initial capillary pressure
1476  feVolInitialPc->getOpPtrVector().push_back(
1477  new OpEvaluateInitiallHead(*this, "VALUES"));
1478 
1479  // set residual face from Neumann terms, i.e. applied pressure
1480  feFaceRhs->getOpPtrVector().push_back(
1481  new OpRhsBcOnValues(*this, "FLUXES", scaleMethodValue));
1482  // set residual finite element operators
1483  headRateAtGaussPts = boost::make_shared<VectorDouble>();
1484  // resAtGaussPts = boost::make_shared<VectorDouble>();
1485  feVolRhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1486  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1487  feVolRhs->getOpPtrVector().push_back(
1488  new OpValuesAtGaussPts(*this, "VALUES"));
1489  feVolRhs->getOpPtrVector().push_back(
1490  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1491  feVolRhs->getOpPtrVector().push_back(new OpResidualFlux(*this, "FLUXES"));
1492  feVolRhs->getOpPtrVector().push_back(new OpResidualMass(*this, "VALUES"));
1493  feVolRhs->getOpPtrVector().back().opType =
1494  ForcesAndSourcesCore::UserDataOperator::OPROW;
1495  // set tangent matrix finite element operators
1496  feVolLhs->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
1497  string("VALUES") + "_t", headRateAtGaussPts, MBTET));
1498  // feVolLhs->getOpPtrVector().push_back(
1499  // new
1500  // OpCalculateScalarFieldValues(string("FLUXES")+"_residual",resAtGaussPts,MBTET)
1501  // );
1502  feVolLhs->getOpPtrVector().push_back(
1503  new OpValuesAtGaussPts(*this, "VALUES"));
1504  feVolLhs->getOpPtrVector().push_back(
1505  new OpFluxDivergenceAtGaussPts(*this, "FLUXES"));
1506  feVolLhs->getOpPtrVector().push_back(
1507  new OpTauDotSigma_HdivHdiv(*this, "FLUXES"));
1508  feVolLhs->getOpPtrVector().push_back(new OpVU_L2L2(*this, "VALUES"));
1509  feVolLhs->getOpPtrVector().push_back(
1510  new OpVDivSigma_L2Hdiv(*this, "VALUES", "FLUXES"));
1511  feVolLhs->getOpPtrVector().push_back(
1512  new OpDivTauU_HdivL2(*this, "FLUXES", "VALUES"));
1513 
1514  // Adding finite elements to DM, time solver will ask for it to assemble
1515  // tangent matrices and residuals
1516  boost::shared_ptr<FEMethod> null;
1517  CHKERR DMMoFEMTSSetIFunction(dM, "MIX_BCVALUE", feFaceRhs, null, null);
1518  CHKERR DMMoFEMTSSetIFunction(dM, "MIX", feVolRhs, null, null);
1519  CHKERR DMMoFEMTSSetIJacobian(dM, "MIX", feVolLhs, null, null);
1520 
1521  // setting up post-processing
1522  boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_process(
1524  CHKERR post_process->generateReferenceElementMesh();
1525  CHKERR post_process->addFieldValuesPostProc("VALUES");
1526  CHKERR post_process->addFieldValuesPostProc("VALUES_t");
1527  // CHKERR post_process->addFieldValuesPostProc("FLUXES_residual");
1528  CHKERR post_process->addFieldValuesPostProc("FLUXES");
1529  post_process->getOpPtrVector().push_back(
1530  new OpValuesAtGaussPts(*this, "VALUES"));
1531  post_process->getOpPtrVector().push_back(
1532  new OpPostProcMaterial(*this, post_process->postProcMesh,
1533  post_process->mapGaussPts, "VALUES"));
1534 
1535  // Integrate fluxes on boundary
1536  boost::shared_ptr<ForcesAndSourcesCore> flux_integrate;
1537  flux_integrate = boost::shared_ptr<ForcesAndSourcesCore>(
1539  flux_integrate->getOpPtrVector().push_back(
1540  new OpIntegrateFluxes(*this, "FLUXES"));
1541  int frequency = 1;
1542  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Monitor post proc", "none");
1543  CHKERR PetscOptionsInt(
1544  "-how_often_output",
1545  "frequency how often results are dumped on hard disk", "", frequency,
1546  &frequency, NULL);
1547  ierr = PetscOptionsEnd();
1548  CHKERRG(ierr);
1549 
1550  tsMonitor = boost::shared_ptr<FEMethod>(
1551  new MonitorPostProc(*this, post_process, flux_integrate, frequency));
1552  TsCtx *ts_ctx;
1553  DMMoFEMGetTsCtx(dM, &ts_ctx);
1554  ts_ctx->get_postProcess_to_do_Monitor().push_back(tsMonitor);
1556  }
1557 
1558  Vec D1; ///< Vector with inital head capillary pressure
1559  Vec ghostFlux; ///< Ghost Vector of integrated fluxes
1560 
1561  /**
1562  * \brief Create vectors and matrices
1563  * @return Error code
1564  */
1567  CHKERR DMCreateMatrix(dM, &Aij);
1568  CHKERR DMCreateGlobalVector(dM, &D0);
1569  CHKERR VecDuplicate(D0, &D1);
1570  CHKERR VecDuplicate(D0, &D);
1571  CHKERR VecDuplicate(D0, &F);
1572  int ghosts[] = {0};
1573  int nb_locals = mField.get_comm_rank() == 0 ? 1 : 0;
1574  int nb_ghosts = mField.get_comm_rank() > 0 ? 1 : 0;
1575  CHKERR VecCreateGhost(PETSC_COMM_WORLD, nb_locals, 1, nb_ghosts, ghosts,
1576  &ghostFlux);
1578  }
1579 
1580  /**
1581  * \brief Delete matrices and vector when no longer needed
1582  * @return error code
1583  */
1586  CHKERR MatDestroy(&Aij);
1587  CHKERR VecDestroy(&D);
1588  CHKERR VecDestroy(&D0);
1589  CHKERR VecDestroy(&D1);
1590  CHKERR VecDestroy(&F);
1591  CHKERR VecDestroy(&ghostFlux);
1593  }
1594 
1595  /**
1596  * \brief Calculate boundary conditions for fluxes
1597  * @return Error code
1598  */
1601  // clear vectors
1602  CHKERR VecZeroEntries(D0);
1603  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_FORWARD);
1604  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_FORWARD);
1605  // clear essential bc indices, it could have dofs from other mesh refinement
1606  bcIndices.clear();
1607  // set operator to calculate essential boundary conditions
1608  CHKERR DMoFEMLoopFiniteElements(dM, "MIX_BCFLUX", feFaceBc);
1609  CHKERR VecGhostUpdateBegin(D0, INSERT_VALUES, SCATTER_REVERSE);
1610  CHKERR VecGhostUpdateEnd(D0, INSERT_VALUES, SCATTER_REVERSE);
1611  CHKERR VecAssemblyBegin(D0);
1612  CHKERR VecAssemblyEnd(D0);
1613  double norm2D0;
1614  CHKERR VecNorm(D0, NORM_2, &norm2D0);
1615  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1616  PetscPrintf(PETSC_COMM_WORLD, "norm2D0 = %6.4e\n", norm2D0);
1618  }
1619 
1620  /**
1621  * \brief Calculate inital pressure head distribution
1622  * @return Error code
1623  */
1626  // clear vectors
1627  CHKERR VecZeroEntries(D1);
1628  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_FORWARD);
1629  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_FORWARD);
1630  // Calculate initial pressure head on each element
1632  // Assemble vector
1633  CHKERR VecGhostUpdateBegin(D1, INSERT_VALUES, SCATTER_REVERSE);
1634  CHKERR VecGhostUpdateEnd(D1, INSERT_VALUES, SCATTER_REVERSE);
1635  CHKERR VecAssemblyBegin(D1);
1636  CHKERR VecAssemblyEnd(D1);
1637  // Calculate norm
1638  double norm2D1;
1639  CHKERR VecNorm(D1, NORM_2, &norm2D1);
1640  // CHKERR VecView(D0,PETSC_VIEWER_STDOUT_WORLD);
1641  PetscPrintf(PETSC_COMM_WORLD, "norm2D1 = %6.4e\n", norm2D1);
1643  }
1644 
1645  /**
1646  * \brief solve problem
1647  * @return error code
1648  */
1649  MoFEMErrorCode solveProblem(bool set_initial_pc = true) {
1651  if (set_initial_pc) {
1652  // Set initial head
1653  CHKERR DMoFEMMeshToLocalVector(dM, D1, INSERT_VALUES, SCATTER_REVERSE);
1654  }
1655 
1656  // Initiate vector from data on the mesh
1657  CHKERR DMoFEMMeshToLocalVector(dM, D, INSERT_VALUES, SCATTER_FORWARD);
1658 
1659  // Create time solver
1660  TS ts;
1661  CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
1662  // Use backward Euler method
1663  CHKERR TSSetType(ts, TSBEULER);
1664  // Set final time
1665  double ftime = 1;
1666  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
1667  // Setup solver from commabd line
1668  CHKERR TSSetFromOptions(ts);
1669  // Set DM to TS
1670  CHKERR TSSetDM(ts, dM);
1671 #if PETSC_VERSION_GE(3, 7, 0)
1672  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1673 #endif
1674  // Set-up monitor
1675  TsCtx *ts_ctx;
1676  DMMoFEMGetTsCtx(dM, &ts_ctx);
1677  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx, PETSC_NULL);
1678 
1679  // This add SNES monitor, to show error by fields. It is dirty trick
1680  // to add monitor, so code is hidden from doxygen
1681  CHKERR TSSetSolution(ts, D);
1682  CHKERR TSSetUp(ts);
1683  SNES snes;
1684  CHKERR TSGetSNES(ts, &snes);
1685 
1686 #if PETSC_VERSION_GE(3, 7, 0)
1687  {
1688  PetscViewerAndFormat *vf;
1689  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1690  PETSC_VIEWER_DEFAULT, &vf);
1691  CHKERR SNESMonitorSet(
1692  snes,
1693  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1694  void *))SNESMonitorFields,
1695  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1696  }
1697 #else
1698  {
1699  CHKERR SNESMonitorSet(snes,
1700  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
1701  void *))SNESMonitorFields,
1702  0, 0);
1703  }
1704 #endif
1705 
1706  CHKERR TSSolve(ts, D);
1707 
1708  // Get statistic form TS and print it
1709  CHKERR TSGetTime(ts, &ftime);
1710  PetscInt steps, snesfails, rejects, nonlinits, linits;
1711  CHKERR TSGetTimeStepNumber(ts, &steps);
1712  CHKERR TSGetSNESFailures(ts, &snesfails);
1713  CHKERR TSGetStepRejections(ts, &rejects);
1714  CHKERR TSGetSNESIterations(ts, &nonlinits);
1715  CHKERR TSGetKSPIterations(ts, &linits);
1716  PetscPrintf(PETSC_COMM_WORLD,
1717  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
1718  "%D, linits %D\n",
1719  steps, rejects, snesfails, ftime, nonlinits, linits);
1720 
1722  }
1723 };
1724 
1725 } // namespace MixTransport
1726 
1727 #endif // __UNSATURATD_FLOW_HPP__
static double scaleZ
Scale z direction.
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
Deprecated interface functions.
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:179
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:445
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:425
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
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.
Mat & ts_B
Preconditioner for ts_A.
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:483
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:550
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:202
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:514
OpResidualFlux(UnsaturatedFlowElement &ctx, const std::string &flux_name)
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:772
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:507
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
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)
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:67
virtual MoFEMErrorCode calSe()
auto getFTensor0IntegrationWeight()
Get integration weights.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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:105
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:602
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.
constexpr int order
Vec & ts_F
residual vector
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
virtual double initialPcEval() const =0
Initialize head.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
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:1001
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.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
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:719
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:413
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:73
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.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
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:982
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.
MoFEMErrorCode buildProblem(Range zero_flux_ents, BitRefLevel ref_level=BitRefLevel().set(0))
Build problem.
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
double w(const double g, const double t)
Definition: ContactOps.hpp:160
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:180