v0.10.0
AnalyticalDirichlet.hpp
Go to the documentation of this file.
1 /** \file AnalyticalDirichlet.hpp
2 
3  Enforce Dirichlet boundary condition for given analytical function,
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifndef __ANALYTICALDIRICHLETBC_HPP__
22 #define __ANALYTICALDIRICHLETBC_HPP__
23 
24 using namespace boost::numeric;
25 using namespace MoFEM;
26 
27 /** \brief Analytical Dirichlet boundary conditions
28  \ingroup user_modules
29  */
31 
32  /** \brief finite element to approximate analytical solution on surface
33  */
34  struct ApproxField {
35 
37 
38  int addToRule; ///< this is add to integration rule if 2nd order geometry
39  ///< approximation
41  : MoFEM::FaceElementForcesAndSourcesCore(m_field), addToRule(4) {}
42  int getRule(int order) { return 2 * order + addToRule; };
43  };
44 
45  ApproxField(MoFEM::Interface &m_field) : feApprox(m_field) {}
46  virtual ~ApproxField() {}
47 
49  MyTriFE &getLoopFeApprox() { return feApprox; }
50 
52  struct OpHoCoord
54 
56  OpHoCoord(const std::string field_name, MatrixDouble &ho_coords);
57 
58  MoFEMErrorCode doWork(int side, EntityType type,
60  };
61 
62  /** \brief Lhs operator used to build matrix
63  */
64  struct OpLhs
66 
68  OpLhs(const std::string field_name, MatrixDouble &ho_coords);
69 
71  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
72  EntityType col_type,
75  };
76 
77  /** \brief Rhs operator used to build matrix
78  */
79  template <typename FUNEVAL>
80  struct OpRhs
82 
83  // Range tRis;
85  boost::shared_ptr<FUNEVAL> functionEvaluator;
87 
88  OpRhs(const std::string field_name,
89  // Range tris,
90  MatrixDouble &ho_coords,
91  boost::shared_ptr<FUNEVAL> function_evaluator, int field_number)
93  field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
94  hoCoords(ho_coords), functionEvaluator(function_evaluator),
95  fieldNumber(field_number) {}
96 
99 
100  MoFEMErrorCode doWork(int side, EntityType type,
103 
104  unsigned int nb_row = data.getIndices().size();
105  if (nb_row == 0)
107 
108  const auto &dof_ptr = data.getFieldDofs()[0];
109  unsigned int rank = dof_ptr->getNbOfCoeffs();
110 
111  NTf.resize(nb_row / rank);
112  iNdices.resize(nb_row / rank);
113 
114  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
115 
116  double x, y, z;
117  double val = getGaussPts()(2, gg);
118  if (hoCoords.size1() == data.getN().size1()) {
119  double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
120  val *= area;
121  x = hoCoords(gg, 0);
122  y = hoCoords(gg, 1);
123  z = hoCoords(gg, 2);
124  } else {
125  val *= getArea();
126  x = getCoordsAtGaussPts()(gg, 0);
127  y = getCoordsAtGaussPts()(gg, 1);
128  z = getCoordsAtGaussPts()(gg, 2);
129  }
130 
131  VectorDouble a = (*functionEvaluator)(x, y, z)[fieldNumber];
132 
133  if (a.size() != rank) {
134  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
135  "data inconsistency");
136  }
137 
138  for (unsigned int rr = 0; rr < rank; rr++) {
139 
140  ublas::noalias(iNdices) = ublas::vector_slice<VectorInt>(
141  data.getIndices(),
142  ublas::slice(rr, rank, data.getIndices().size() / rank));
143 
144  noalias(NTf) = data.getN(gg, nb_row / rank) * a[rr] * val;
145  ierr = VecSetValues(getFEMethod()->snes_f, iNdices.size(),
146  &iNdices[0], &*NTf.data().begin(), ADD_VALUES);
147  CHKERRG(ierr);
148  }
149  }
150 
152  }
153  };
154  };
155 
156  /**
157  * \brief Structure used to enforce analytical boundary conditions
158  */
160 
161  DirichletBC(MoFEM::Interface &m_field, const std::string &field, Mat A,
162  Vec X, Vec F);
163 
164  DirichletBC(MoFEM::Interface &m_field, const std::string &field);
165 
166  boost::shared_ptr<Range> trisPtr;
167 
168  MoFEMErrorCode iNitalize();
169  MoFEMErrorCode iNitalize(Range &tris);
170  };
171 
174 
175  /**
176  * \brief Set operators used to calculate the rhs vector and the lhs matrix
177 
178  * To enforce analytical function on boundary, first function has to be
179  approximated
180  * by finite element base functions. This is done by solving system of linear
181  * equations, Following function set finite element operators to calculate
182  * the left hand side matrix and the left hand side matrix.
183 
184  * @param m_field interface
185  * @param field_name field name
186  * @param function_evaluator analytical function to evaluate
187  * @param field_number field index
188  * @param nodals_positions name of the field for ho-geometry description
189  * @return error code
190  */
191  template <typename FUNEVAL>
193  setApproxOps(MoFEM::Interface &m_field, const std::string field_name,
194  boost::shared_ptr<FUNEVAL> function_evaluator,
195  const int field_number = 0,
196  const string nodals_positions = "MESH_NODE_POSITIONS") {
198  if (approxField.getLoopFeApprox().getOpPtrVector().empty()) {
199  if (m_field.check_field(nodals_positions)) {
200  approxField.getLoopFeApprox().getOpPtrVector().push_back(
201  new ApproxField::OpHoCoord(nodals_positions, approxField.hoCoords));
202  }
203  approxField.getLoopFeApprox().getOpPtrVector().push_back(
204  new ApproxField::OpLhs(field_name, approxField.hoCoords));
205  }
206  approxField.getLoopFeApprox().getOpPtrVector().push_back(
207  new ApproxField::OpRhs<FUNEVAL>(field_name, approxField.hoCoords,
208  function_evaluator, field_number));
210  }
211 
212  // /**
213  // * \deprecated no need to use function with argument of triangle range
214  // */
215  // template<typename FUNEVAL> DEPRECATED MoFEMErrorCode setApproxOps(
216  // MoFEM::Interface &m_field,
217  // string field_name,
218  // Range& tris,
219  // boost::shared_ptr<FUNEVAL> function_evaluator,
220  // int field_number = 0,
221  // string nodals_positions = "MESH_NODE_POSITIONS"
222  // ) {
223  // return setApproxOps(
224  // m_field,field_name,tris,function_evaluator,field_number,nodals_positions
225  // );
226  // }
227 
228  /**
229  * \brief set finite element
230  * @param m_field mofem interface
231  * @param fe finite element name
232  * @param field field name
233  * @param tris faces where analytical boundary is given
234  * @param nodals_positions field having higher order geometry description
235  * @return error code
236  */
238  setFiniteElement(MoFEM::Interface &m_field, string fe, string field,
239  Range &tris,
240  string nodals_positions = "MESH_NODE_POSITIONS");
241 
242  // /**
243  // \deprecated use setFiniteElement instead
244  // */
245  // DEPRECATED MoFEMErrorCode initializeProblem(
246  // MoFEM::Interface &m_field,
247  // string fe,
248  // string field,
249  // Range& tris,
250  // string nodals_positions = "MESH_NODE_POSITIONS"
251  // ) {
252  // return setFiniteElement(m_field,fe,field,tris,nodals_positions);
253  // }
254 
255  Mat A;
256  Vec D, F;
258 
259  /**
260  * \brief set problem solver and create matrices and vectors
261  * @param m_field mofem interface
262  * @param problem problem name
263  * @return error code
264  */
265  MoFEMErrorCode setUpProblem(MoFEM::Interface &m_field, string problem);
266 
267  // /**
268  // * \deprecated use setUpProblem instead
269  // */
270  // DEPRECATED MoFEMErrorCode setProblem(
271  // MoFEM::Interface &m_field,string problem
272  // ) {
273  // return setUpProblem(m_field,problem);
274  // }
275 
276  /**
277  * \brief solve boundary problem
278  *
279 
280  * This functions solve for DOFs on boundary where analytical solution is
281  * given. i.e. finding values of DOFs which approximate analytical solution.
282 
283  * @param m_field mofem interface
284  * @param problem problem name
285  * @param fe finite element name
286  * @param bc Driblet boundary structure used to apply boundary
287  conditions
288  * @param tris triangles on boundary
289  * @return error code
290  */
291  MoFEMErrorCode solveProblem(MoFEM::Interface &m_field, string problem,
292  string fe, DirichletBC &bc, Range &tris);
293 
294  /**
295  * \brief solve boundary problem
296  *
297 
298  * This functions solve for DOFs on boundary where analytical solution is
299  * given.
300 
301  * @param m_field mofem interface
302  * @param problem problem name
303  * @param fe finite element name
304  * @param bc Driblet boundary structure used to apply boundary
305  conditions
306  * @return [description]
307  */
308  MoFEMErrorCode solveProblem(MoFEM::Interface &m_field, string problem,
309  string fe, DirichletBC &bc);
310 
311  /**
312  * \brief Destroy problem
313  *
314  * Destroy matrices and vectors used to solve boundary problem, i.e. finding
315  * values of DOFs which approximate analytical solution.
316  *
317  * @return error code
318  */
319  MoFEMErrorCode destroyProblem();
320 };
321 
322 #endif //__ANALYTICALDIRICHLETBC_HPP__
Deprecated interface functions.
ApproxField(MoFEM::Interface &m_field)
const double D
diffusivity
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
Structure used to enforce analytical boundary conditions.
Rhs operator used to build matrix.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
Analytical Dirichlet boundary conditions.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
const VectorInt & getIndices() const
Get global indices of dofs on entity.
OpRhs(const std::string field_name, MatrixDouble &ho_coords, boost::shared_ptr< FUNEVAL > function_evaluator, int field_number)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
MoFEMErrorCode setApproxOps(MoFEM::Interface &m_field, const std::string field_name, boost::shared_ptr< FUNEVAL > function_evaluator, const int field_number=0, const string nodals_positions="MESH_NODE_POSITIONS")
Set operators used to calculate the rhs vector and the lhs matrix.
constexpr int order
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
virtual bool check_field(const std::string &name) const =0
check if field is in database
Data on single entity (This is passed as argument to DataOperator::doWork)
structure to get information form mofem into DataForcesAndSourcesCore
boost::shared_ptr< FUNEVAL > functionEvaluator
finite element to approximate analytical solution on surface
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
Lhs operator used to build matrix.