v0.14.0
Loading...
Searching...
No Matches
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
8
9#ifndef __ANALYTICALDIRICHLETBC_HPP__
10#define __ANALYTICALDIRICHLETBC_HPP__
11
12using namespace boost::numeric;
13using namespace MoFEM;
14
15/** \brief Analytical Dirichlet boundary conditions
16 \ingroup user_modules
17 */
19
20 /** \brief finite element to approximate analytical solution on surface
21 */
22 struct ApproxField {
23
25
26 int addToRule; ///< this is add to integration rule if 2nd order geometry
27 ///< approximation
30 int getRule(int order) { return 2 * order + addToRule; };
31 };
32
33 ApproxField(MoFEM::Interface &m_field) : feApprox(m_field) {}
34 virtual ~ApproxField() = default;
35
38
39 /** \brief Lhs operator used to build matrix
40 */
41 struct OpLhs
43
44 OpLhs(const std::string field_name);
45
47 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
48 EntityType col_type,
51 };
52
53 /** \brief Rhs operator used to build matrix
54 */
55 template <typename FUNEVAL>
56 struct OpRhs
58
59 // Range tRis;
60 boost::shared_ptr<FUNEVAL> functionEvaluator;
62
63 OpRhs(const std::string field_name,
64 boost::shared_ptr<FUNEVAL> function_evaluator, int field_number)
67 functionEvaluator(function_evaluator), fieldNumber(field_number) {}
68
71
75
76 unsigned int nb_row = data.getIndices().size();
77 if (nb_row == 0)
79
80 const auto &dof_ptr = data.getFieldDofs()[0];
81 unsigned int rank = dof_ptr->getNbOfCoeffs();
82
83 const auto &gauss_pts = getGaussPts();
84 const auto &coords_at_gauss_pts = getCoordsAtGaussPts();
85
86 NTf.resize(nb_row / rank);
87 iNdices.resize(nb_row / rank);
88
89 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
90
91 const double area = norm_2(getNormalsAtGaussPts(gg)) * 0.5;
92 const double val = gauss_pts(2, gg) * area;
93 const double x = coords_at_gauss_pts(gg, 0);
94 const double y = coords_at_gauss_pts(gg, 1);
95 const double z = coords_at_gauss_pts(gg, 2);
96
97 VectorDouble a = (*functionEvaluator)(x, y, z)[fieldNumber];
98
99 if (a.size() != rank) {
100 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
101 "data inconsistency");
102 }
103
104 for (unsigned int rr = 0; rr < rank; rr++) {
105
106 ublas::noalias(iNdices) = ublas::vector_slice<VectorInt>(
107 data.getIndices(),
108 ublas::slice(rr, rank, data.getIndices().size() / rank));
109
110 noalias(NTf) = data.getN(gg, nb_row / rank) * a[rr] * val;
111 CHKERR VecSetValues(getFEMethod()->snes_f, iNdices.size(),
112 &iNdices[0], &*NTf.data().begin(), ADD_VALUES);
113 }
114 }
115
117 }
118 };
119 };
120
121 /**
122 * \brief Structure used to enforce analytical boundary conditions
123 */
125
126 DirichletBC(MoFEM::Interface &m_field, const std::string &field, Mat A,
127 Vec X, Vec F);
128
129 DirichletBC(MoFEM::Interface &m_field, const std::string &field);
130
131 boost::shared_ptr<Range> trisPtr;
132
135 };
136
139
140 /**
141 * \brief Set operators used to calculate the rhs vector and the lhs matrix
142
143 * To enforce analytical function on boundary, first function has to be
144 approximated
145 * by finite element base functions. This is done by solving system of linear
146 * equations, Following function set finite element operators to calculate
147 * the left hand side matrix and the left hand side matrix.
148
149 * @param m_field interface
150 * @param field_name field name
151 * @param function_evaluator analytical function to evaluate
152 * @param field_number field index
153 * @param nodals_positions name of the field for ho-geometry description
154 * @return error code
155 */
156 template <typename FUNEVAL>
158 setApproxOps(MoFEM::Interface &m_field, const std::string field_name,
159 boost::shared_ptr<FUNEVAL> function_evaluator,
160 const int field_number = 0,
161 const string nodals_positions = "MESH_NODE_POSITIONS") {
164 if (m_field.check_field(nodals_positions))
166 false, false);
169 }
171 new ApproxField::OpRhs<FUNEVAL>(field_name, function_evaluator,
172 field_number));
174 }
175
176 // /**
177 // * \deprecated no need to use function with argument of triangle range
178 // */
179 // template<typename FUNEVAL> DEPRECATED MoFEMErrorCode setApproxOps(
180 // MoFEM::Interface &m_field,
181 // string field_name,
182 // Range& tris,
183 // boost::shared_ptr<FUNEVAL> function_evaluator,
184 // int field_number = 0,
185 // string nodals_positions = "MESH_NODE_POSITIONS"
186 // ) {
187 // return setApproxOps(
188 // m_field,field_name,tris,function_evaluator,field_number,nodals_positions
189 // );
190 // }
191
192 /**
193 * \brief set finite element
194 * @param m_field mofem interface
195 * @param fe finite element name
196 * @param field field name
197 * @param tris faces where analytical boundary is given
198 * @param nodals_positions field having higher order geometry description
199 * @return error code
200 */
202 setFiniteElement(MoFEM::Interface &m_field, string fe, string field,
203 Range &tris,
204 string nodals_positions = "MESH_NODE_POSITIONS");
205
206 // /**
207 // \deprecated use setFiniteElement instead
208 // */
209 // DEPRECATED MoFEMErrorCode initializeProblem(
210 // MoFEM::Interface &m_field,
211 // string fe,
212 // string field,
213 // Range& tris,
214 // string nodals_positions = "MESH_NODE_POSITIONS"
215 // ) {
216 // return setFiniteElement(m_field,fe,field,tris,nodals_positions);
217 // }
218
219 Mat A;
220 Vec D, F;
222
223 /**
224 * \brief set problem solver and create matrices and vectors
225 * @param m_field mofem interface
226 * @param problem problem name
227 * @return error code
228 */
229 MoFEMErrorCode setUpProblem(MoFEM::Interface &m_field, string problem);
230
231 // /**
232 // * \deprecated use setUpProblem instead
233 // */
234 // DEPRECATED MoFEMErrorCode setProblem(
235 // MoFEM::Interface &m_field,string problem
236 // ) {
237 // return setUpProblem(m_field,problem);
238 // }
239
240 /**
241 * \brief solve boundary problem
242 *
243
244 * This functions solve for DOFs on boundary where analytical solution is
245 * given. i.e. finding values of DOFs which approximate analytical solution.
246
247 * @param m_field mofem interface
248 * @param problem problem name
249 * @param fe finite element name
250 * @param bc Driblet boundary structure used to apply boundary
251 conditions
252 * @param tris triangles on boundary
253 * @return error code
254 */
255 MoFEMErrorCode solveProblem(MoFEM::Interface &m_field, string problem,
256 string fe, DirichletBC &bc, Range &tris);
257
258 /**
259 * \brief solve boundary problem
260 *
261
262 * This functions solve for DOFs on boundary where analytical solution is
263 * given.
264
265 * @param m_field mofem interface
266 * @param problem problem name
267 * @param fe finite element name
268 * @param bc Driblet boundary structure used to apply boundary
269 conditions
270 * @return [description]
271 */
272 MoFEMErrorCode solveProblem(MoFEM::Interface &m_field, string problem,
273 string fe, DirichletBC &bc);
274
275 /**
276 * \brief Destroy problem
277 *
278 * Destroy matrices and vectors used to solve boundary problem, i.e. finding
279 * values of DOFs which approximate analytical solution.
280 *
281 * @return error code
282 */
284};
285
286#endif //__ANALYTICALDIRICHLETBC_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
virtual bool check_field(const std::string &name) const =0
check if field is in database
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr auto field_name
Lhs operator used to build matrix.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Rhs operator used to build matrix.
OpRhs(const std::string field_name, boost::shared_ptr< FUNEVAL > function_evaluator, int field_number)
boost::shared_ptr< FUNEVAL > functionEvaluator
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
finite element to approximate analytical solution on surface
ApproxField(MoFEM::Interface &m_field)
Structure used to enforce analytical boundary conditions.
Analytical Dirichlet boundary conditions.
MoFEMErrorCode destroyProblem()
Destroy problem.
MoFEMErrorCode setUpProblem(MoFEM::Interface &m_field, string problem)
set problem solver and create matrices and vectors
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.
MoFEMErrorCode setFiniteElement(MoFEM::Interface &m_field, string fe, string field, Range &tris, string nodals_positions="MESH_NODE_POSITIONS")
set finite element
MoFEMErrorCode solveProblem(MoFEM::Interface &m_field, string problem, string fe, DirichletBC &bc, Range &tris)
solve boundary problem
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.