v0.15.0
Loading...
Searching...
No Matches
EssentialDisplacementCubitBcData.hpp
Go to the documentation of this file.
1/**
2 * @file EssentialDisplacementCubitBcData.hpp
3 * @brief Implementations specialisation for DisplacementCubitBcData
4 * @version 0.13.2
5 * @date 2022-09-18
6 *
7 * @copyright Copyright (c) 2022
8 *
9 */
10
11#ifndef _ESSENTIAL_DISPLACEMENTCUBITBCDATA_HPP__
12#define _ESSENTIAL_DISPLACEMENTCUBITBCDATA_HPP__
13
14namespace MoFEM {
15
16/**
17 * @brief A specialized version of DisplacementCubitBcData that includes an
18 * additional rotation offset.
19 *
20 * DisplacementCubitBcDataWithRotation extends the DisplacementCubitBcData
21 * structure with a rotation offset and a method to calculate the rotated
22 * displacement given the rotation angles and coordinates.
23 */
25
26 std::array<double, 3> rotOffset;
28
29 /**
30 * @brief Calculates the rotated displacement given the rotation angles,
31 * coordinates, and an optional offset.
32 *
33 * @param angles A FTensor::Tensor1 containing the rotation angles.
34 * @param coordinates A FTensor::Tensor1 containing the coordinates.
35 * @param offset An optional FTensor::Tensor1 containing the offset
36 * (default is {0., 0., 0.}).
37 * @return FTensor::Tensor1<double, 3> representing the rotated displacement.
38 */
43 0., 0., 0.}) {
44
45 FTensor::Index<'i', 3> i;
46 FTensor::Index<'j', 3> j;
47 FTensor::Index<'k', 3> k;
48
49 FTensor::Tensor1<double, 3> rotated_displacement;
50 rotated_displacement(i) = 0;
51
52 auto get_rotation = [&](auto &omega) {
53 FTensor::Tensor2<double, 3, 3> rotation_matrix;
54
55 constexpr auto kronecker_delta = FTensor::Kronecker_Delta<int>();
56 rotation_matrix(i, j) = kronecker_delta(i, j);
57
58 const double angle = sqrt(omega(i) * omega(i));
59 if (std::abs(angle) < std::numeric_limits<double>::epsilon())
60 return rotation_matrix;
61
63 t_omega(i, j) = FTensor::levi_civita<double>(i, j, k) * omega(k);
64 const double a = sin(angle) / angle;
65 const double sin_squared = sin(angle / 2.) * sin(angle / 2.);
66 const double b = 2. * sin_squared / (angle * angle);
67 rotation_matrix(i, j) += a * t_omega(i, j);
68 rotation_matrix(i, j) += b * t_omega(i, k) * t_omega(k, j);
69
70 return rotation_matrix;
71 };
72
73 auto rotation = get_rotation(angles);
74 FTensor::Tensor1<double, 3> coordinate_distance;
75 coordinate_distance(i) = offset(i) - coordinates(i);
76 rotated_displacement(i) =
77 coordinate_distance(i) - rotation(j, i) * coordinate_distance(j);
78
79 return rotated_displacement;
80 }
81};
82
83/**
84 * @brief Specialization for DisplacementCubitBcData
85 *
86 * Specialization to enforce blocksets which DisplacementCubitBcData ptr. That
87 * is to enforce displacement constraints. set
88 *
89 * @tparam
90 */
93 boost::shared_ptr<FEMethod> fe_ptr,
94 std::vector<boost::shared_ptr<ScalingMethod>> smv,
95 bool get_coords = false);
96
97 virtual ~EssentialPreProc() = default;
98
99 MoFEMErrorCode operator()();
100
101protected:
103 boost::weak_ptr<FEMethod> fePtr;
106};
107
108/**
109 * @brief Specialization for DisplacementCubitBcData
110 *
111 * @tparam
112 */
115 boost::shared_ptr<FEMethod> fe_ptr, double diag,
116 SmartPetscObj<Vec> rhs = nullptr);
117
118 virtual ~EssentialPostProcRhs() = default;
119
120 MoFEMErrorCode operator()();
121
122protected:
124 boost::weak_ptr<FEMethod> fePtr;
125 double vDiag;
127};
128
129/**
130 * @brief Specialization for DisplacementCubitBcData
131 *
132 * @tparam
133 */
136 boost::shared_ptr<FEMethod> fe_ptr, double diag,
137 SmartPetscObj<Mat> lhs = nullptr,
138 SmartPetscObj<AO> ao = nullptr);
139
140 virtual ~EssentialPostProcLhs() = default;
141
142 MoFEMErrorCode operator()();
143
144protected:
146 boost::weak_ptr<FEMethod> fePtr;
147 double vDiag;
150};
151
152/**
153 * @brief Specialization for DisplacementCubitBcData
154 *
155 * @tparam
156 */
159 MoFEM::Interface &m_field, boost::shared_ptr<FEMethod> fe_ptr,
160 SmartPetscObj<Vec> rhs = nullptr,
161 LogManager::SeverityLevel sev = Sev::inform,
162 boost::shared_ptr<std::vector<double>> reaction_ptr = nullptr);
163
164 MoFEMErrorCode operator()();
165
166protected:
168 boost::weak_ptr<FEMethod> fePtr;
171 PetscBool printBlockName; //< print block name when printing reaction
172 boost::shared_ptr<std::vector<double>> reactionPtr;
173 char reactionBlockName[255];
174};
175
176template <int FIELD_DIM, AssemblyType A, IntegrationType I, typename OpBase>
178 : FormsIntegrators<OpBase>::template Assembly<A>::template LinearForm<
179 I>::template OpSource<1, FIELD_DIM> {
180
183
184 OpEssentialRhsImpl(const std::string field_name,
185 boost::shared_ptr<DisplacementCubitBcData> bc_data,
186 boost::shared_ptr<Range> ents_ptr,
187 std::vector<boost::shared_ptr<ScalingMethod>> smv);
188
189private:
194
196 double z) {
199 tAngles, FTensor::Tensor1<double, 3>{x, y, z}, tOffset);
201 t_ret(i) = tVal(i) + rot(i);
202 return t_ret;
203 };
204
206 initDispFunction(boost::shared_ptr<DisplacementCubitBcData> bc_data) {
207 if (bc_data->data.flag4 || bc_data->data.flag5 || bc_data->data.flag6) {
208 return [this](double x, double y, double z) {
209 return this->rotFunction(x, y, z);
210 };
211 }
212 return [this](double x, double y, double z) { return tVal; };
213 }
214};
215
216template <int FIELD_DIM, AssemblyType A, IntegrationType I, typename OpBase>
218 OpEssentialRhsImpl(const std::string field_name,
219 boost::shared_ptr<DisplacementCubitBcData> bc_data,
220 boost::shared_ptr<Range> ents_ptr,
221 std::vector<boost::shared_ptr<ScalingMethod>> smv)
222 : OpSource(field_name, initDispFunction(bc_data), ents_ptr),
223 vecOfTimeScalingMethods(smv) {
224 static_assert(FIELD_DIM > 1, "Is not implemented for scalar field");
225
227 tVal(i) = 0;
228 tAngles(i) = 0;
229 tOffset(i) = 0;
230
231 if (bc_data->data.flag1 == 1)
232 tVal(0) = -bc_data->data.value1;
233 if (bc_data->data.flag2 == 1 && FIELD_DIM > 1)
234 tVal(1) = -bc_data->data.value2;
235 if (bc_data->data.flag3 == 1 && FIELD_DIM > 2)
236 tVal(2) = -bc_data->data.value3;
237 if (bc_data->data.flag4 == 1 && FIELD_DIM > 2)
238 tAngles(0) = -bc_data->data.value4;
239 if (bc_data->data.flag5 == 1 && FIELD_DIM > 2)
240 tAngles(1) = -bc_data->data.value5;
241 if (bc_data->data.flag6 == 1 && FIELD_DIM > 1)
242 tAngles(2) = -bc_data->data.value6;
243
244 if (auto ext_bc_data =
245 dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
246 bc_data.get())) {
247 for (int a = 0; a != 3; ++a)
248 tOffset(a) = ext_bc_data->rotOffset[a];
249 }
250
251 this->timeScalingFun = [this](const double t) {
252 double s = 1;
253 for (auto &o : vecOfTimeScalingMethods) {
254 s *= o->getScale(t);
255 }
256 return s;
257 };
258}
259
260template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
261 typename OpBase>
263 OpBase>
264 : FormsIntegrators<OpBase>::template Assembly<A>::template BiLinearForm<
265 I>::template OpMass<BASE_DIM, FIELD_DIM> {
266
269
270 OpEssentialLhsImpl(const std::string field_name,
271 boost::shared_ptr<Range> ents_ptr);
272};
273
274template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
275 typename OpBase>
277 OpBase>::OpEssentialLhsImpl(const std::string field_name,
278 boost::shared_ptr<Range>
279 ents_ptr)
280 : OpMass(
282
283 [](double, double, double) constexpr { return 1; },
284
285 ents_ptr) {}
286
287template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
288 typename OpBase>
290
292 OpBase>,
293 A, I, OpBase
294
295 > {
296
298
300
301 MoFEM::Interface &m_field,
302 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
303 const std::string problem_name, std::string field_name,
304 boost::shared_ptr<MatrixDouble> field_mat_ptr,
305 std::vector<boost::shared_ptr<ScalingMethod>> smv
306
307 ) {
309
310 using OP =
311 typename EssentialBC<OpBase>::template Assembly<A>::template LinearForm<
312 I>::template OpEssentialRhs<DisplacementCubitBcData, BASE_DIM,
313 FIELD_DIM>;
314 using OpInternal = typename FormsIntegrators<OpBase>::template Assembly<
315 A>::template LinearForm<I>::template OpBaseTimesVector<BASE_DIM,
316 FIELD_DIM, 1>;
317
318 auto add_op = [&](auto &bcs) {
320 for (auto &m : bcs) {
321 if (auto bc = m.second->dispBcPtr) {
322 auto &bc_id = m.first;
323 auto regex_str =
324 (boost::format("%s_%s_(.*)") % problem_name % field_name).str();
325 if (std::regex_match(bc_id, std::regex(regex_str))) {
326 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "OpEssentialRhs") << *bc;
327 pipeline.push_back(
328 new OpSetBc(field_name, false, m.second->getBcMarkersPtr()));
329 pipeline.push_back(
330 new OP(field_name, bc, m.second->getBcEntsPtr(), smv));
331 pipeline.push_back(new OpInternal(
332 field_name, field_mat_ptr,
333 [](double, double, double) constexpr { return 1.; },
334 m.second->getBcEntsPtr()));
335 pipeline.push_back(new OpUnSetBc(field_name));
336 }
337 }
338 }
339 MOFEM_LOG_CHANNEL("SELF");
341 };
342
343 CHKERR add_op(m_field.getInterface<BcManager>()->getBcMapByBlockName());
344
346 }
347};
348
349template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
350 typename OpBase>
352
354 OpBase>,
355 A, I, OpBase
356
357 > {
358
360
362
363 MoFEM::Interface &m_field,
364 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
365 const std::string problem_name, std::string field_name
366
367 ) {
369
370 using OP = typename EssentialBC<OpBase>::template Assembly<A>::
371 template BiLinearForm<I>::template OpEssentialLhs<
373
374 auto add_op = [&](auto &bcs) {
376 for (auto &m : bcs) {
377 if (auto bc = m.second->dispBcPtr) {
378 auto &bc_id = m.first;
379 auto regex_str =
380 (boost::format("%s_%s_(.*)") % problem_name % field_name).str();
381 if (std::regex_match(bc_id, std::regex(regex_str))) {
382 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "OpEssentialLhs") << *bc;
383 pipeline.push_back(
384 new OpSetBc(field_name, false, m.second->getBcMarkersPtr()));
385 pipeline.push_back(new OP(field_name, m.second->getBcEntsPtr()));
386 pipeline.push_back(new OpUnSetBc(field_name));
387 }
388 }
389 }
390 MOFEM_LOG_CHANNEL("SELF");
392 };
393
394 CHKERR add_op(m_field.getInterface<BcManager>()->getBcMapByBlockName());
395
397 }
398};
399
400} // namespace MoFEM
401
402#endif // _ESSENTIAL_DISPLACEMENTCUBITBCDATA_HPP__
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
constexpr double a
constexpr int FIELD_DIM
Kronecker Delta class.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int BASE_DIM
constexpr double omega
Save field DOFS on vertices/tags.
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
boost::function< FTensor::Tensor1< double, DIM >( const double, const double, const double)> VectorFun
Vector function type.
SeverityLevel
Severity levels.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
std::vector< boost::shared_ptr< ScalingMethod > > VecOfTimeScalingMethods
Vector of time scaling methods.
Definition Natural.hpp:20
constexpr IntegrationType I
constexpr AssemblyType A
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
FTensor::Index< 'm', 3 > m
static MoFEMErrorCode add(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, const std::string problem_name, std::string field_name)
Function (factory) for setting operators for lhs pipeline.
static MoFEMErrorCode add(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, const std::string problem_name, std::string field_name, boost::shared_ptr< MatrixDouble > field_mat_ptr, std::vector< boost::shared_ptr< ScalingMethod > > smv)
Function (factory) for setting operators for rhs pipeline.
Definition Essential.hpp:92
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Deprecated interface functions.
A specialized version of DisplacementCubitBcData that includes an additional rotation offset.
static FTensor::Tensor1< double, 3 > GetRotDisp(const FTensor::Tensor1< double, 3 > &angles, FTensor::Tensor1< double, 3 > coordinates, FTensor::Tensor1< double, 3 > offset=FTensor::Tensor1< double, 3 >{ 0., 0., 0.})
Calculates the rotated displacement given the rotation angles, coordinates, and an optional offset.
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Essential boundary conditions.
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
Class (Function) to calculate residual side diagonal.
Definition Essential.hpp:49
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
typename FormsIntegrators< OpBase >::template Assembly< A >::template BiLinearForm< I >::template OpMass< BASE_DIM, FIELD_DIM > OpMass
Enforce essential constrains on lhs.
Definition Essential.hpp:81
VectorFun< FIELD_DIM > initDispFunction(boost::shared_ptr< DisplacementCubitBcData > bc_data)
typename FormsIntegrators< OpBase >::template Assembly< A >::template LinearForm< I >::template OpSource< 1, FIELD_DIM > OpSource
FTensor::Tensor1< double, FIELD_DIM > rotFunction(double x, double y, double z)
Enforce essential constrains on rhs.
Definition Essential.hpp:65
Set indices on entities on finite element.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.