v0.14.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
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 MoFEMErrorCode operator()();
98
99protected:
101 boost::weak_ptr<FEMethod> fePtr;
104};
105
106/**
107 * @brief Specialization for DisplacementCubitBcData
108 *
109 * @tparam
110 */
113 boost::shared_ptr<FEMethod> fe_ptr, double diag,
114 SmartPetscObj<Vec> rhs = nullptr);
115
116 MoFEMErrorCode operator()();
117
118protected:
120 boost::weak_ptr<FEMethod> fePtr;
121 double vDiag;
123};
124
125/**
126 * @brief Specialization for DisplacementCubitBcData
127 *
128 * @tparam
129 */
132 boost::shared_ptr<FEMethod> fe_ptr, double diag,
133 SmartPetscObj<Mat> lhs = nullptr,
134 SmartPetscObj<AO> ao = nullptr);
135
136 MoFEMErrorCode operator()();
137
138protected:
140 boost::weak_ptr<FEMethod> fePtr;
141 double vDiag;
144};
145
146/**
147 * @brief Specialization for DisplacementCubitBcData
148 *
149 * @tparam
150 */
153 boost::shared_ptr<FEMethod> fe_ptr,
154 SmartPetscObj<Vec> rhs = nullptr,
155 LogManager::SeverityLevel sev = Sev::inform);
156
157 MoFEMErrorCode operator()();
158
159protected:
161 boost::weak_ptr<FEMethod> fePtr;
164};
165
166template <int FIELD_DIM, AssemblyType A, IntegrationType I, typename OpBase>
168 : FormsIntegrators<OpBase>::template Assembly<A>::template LinearForm<
169 I>::template OpSource<1, FIELD_DIM> {
170
173
174 OpEssentialRhsImpl(const std::string field_name,
175 boost::shared_ptr<DisplacementCubitBcData> bc_data,
176 boost::shared_ptr<Range> ents_ptr,
177 std::vector<boost::shared_ptr<ScalingMethod>> smv);
178
179private:
185
187 double z) {
190 tAngles, FTensor::Tensor1<double, 3>{x, y, z}, tOffset);
192 t_ret(i) = tVal(i) + rot(i);
193 return t_ret;
194 };
195};
196
197template <int FIELD_DIM, AssemblyType A, IntegrationType I, typename OpBase>
199 OpEssentialRhsImpl(const std::string field_name,
200 boost::shared_ptr<DisplacementCubitBcData> bc_data,
201 boost::shared_ptr<Range> ents_ptr,
202 std::vector<boost::shared_ptr<ScalingMethod>> smv)
203 : OpSource(
204 field_name, [this](double, double, double) { return tVal; },
205 ents_ptr),
206 vecOfTimeScalingMethods(smv) {
207 static_assert(FIELD_DIM > 1, "Is not implemented for scalar field");
208
210 tVal(i) = 0;
211 tAngles(i) = 0;
212 tOffset(i) = 0;
213
214 if (bc_data->data.flag1 == 1)
215 tVal(0) = -bc_data->data.value1;
216 if (bc_data->data.flag2 == 1 && FIELD_DIM > 1)
217 tVal(1) = -bc_data->data.value2;
218 if (bc_data->data.flag3 == 1 && FIELD_DIM > 2)
219 tVal(2) = -bc_data->data.value3;
220 if (bc_data->data.flag4 == 1 && FIELD_DIM > 2)
221 tAngles(0) = -bc_data->data.value4;
222 if (bc_data->data.flag5 == 1 && FIELD_DIM > 2)
223 tAngles(1) = -bc_data->data.value5;
224 if (bc_data->data.flag6 == 1 && FIELD_DIM > 1)
225 tAngles(2) = -bc_data->data.value6;
226
227 if (auto ext_bc_data =
228 dynamic_cast<DisplacementCubitBcDataWithRotation const *>(
229 bc_data.get())) {
230 for (int a = 0; a != 3; ++a)
231 tOffset(a) = ext_bc_data->rotOffset[a];
232 }
233
234 this->timeScalingFun = [this](const double t) {
235 double s = 1;
236 for (auto &o : vecOfTimeScalingMethods) {
237 s *= o->getScale(t);
238 }
239 return s;
240 };
241 if (bc_data->data.flag4 || bc_data->data.flag5 || bc_data->data.flag6) {
242 this->dispFunction = [this](double x, double y, double z) {
243 return this->rotFunction(x, y, z);
244 };
245 } else {
246 // use default set at construction
247 }
248}
249
250template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
251 typename OpBase>
253 OpBase>
254 : FormsIntegrators<OpBase>::template Assembly<A>::template BiLinearForm<
255 I>::template OpMass<BASE_DIM, FIELD_DIM> {
256
259
260 OpEssentialLhsImpl(const std::string field_name,
261 boost::shared_ptr<Range> ents_ptr);
262};
263
264template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
265 typename OpBase>
267 OpBase>::OpEssentialLhsImpl(const std::string field_name,
268 boost::shared_ptr<Range>
269 ents_ptr)
270 : OpMass(
272
273 [](double, double, double) constexpr { return 1; },
274
275 ents_ptr) {}
276
277template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
278 typename OpBase>
280
282 OpBase>,
283 A, I, OpBase
284
285 > {
286
288
290
291 MoFEM::Interface &m_field,
292 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
293 const std::string problem_name, std::string field_name,
294 boost::shared_ptr<MatrixDouble> field_mat_ptr,
295 std::vector<boost::shared_ptr<ScalingMethod>> smv
296
297 ) {
299
300 using OP =
301 typename EssentialBC<OpBase>::template Assembly<A>::template LinearForm<
302 I>::template OpEssentialRhs<DisplacementCubitBcData, BASE_DIM,
303 FIELD_DIM>;
304 using OpInternal = typename FormsIntegrators<OpBase>::template Assembly<
305 A>::template LinearForm<I>::template OpBaseTimesVector<BASE_DIM,
306 FIELD_DIM, 1>;
307
308 auto add_op = [&](auto &bcs) {
310 for (auto &m : bcs) {
311 if (auto bc = m.second->dispBcPtr) {
312 auto &bc_id = m.first;
313 auto regex_str =
314 (boost::format("%s_%s_(.*)") % problem_name % field_name).str();
315 if (std::regex_match(bc_id, std::regex(regex_str))) {
316 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "OpEssentialRhs") << *bc;
317 pipeline.push_back(
318 new OpSetBc(field_name, false, m.second->getBcMarkersPtr()));
319 pipeline.push_back(
320 new OP(field_name, bc, m.second->getBcEntsPtr(), smv));
321 pipeline.push_back(new OpInternal(
322 field_name, field_mat_ptr,
323 [](double, double, double) constexpr { return 1.; },
324 m.second->getBcEntsPtr()));
325 pipeline.push_back(new OpUnSetBc(field_name));
326 }
327 }
328 }
329 MOFEM_LOG_CHANNEL("SELF");
331 };
332
333 CHKERR add_op(m_field.getInterface<BcManager>()->getBcMapByBlockName());
334
336 }
337};
338
339template <int BASE_DIM, int FIELD_DIM, AssemblyType A, IntegrationType I,
340 typename OpBase>
342
344 OpBase>,
345 A, I, OpBase
346
347 > {
348
350
352
353 MoFEM::Interface &m_field,
354 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
355 const std::string problem_name, std::string field_name
356
357 ) {
359
360 using OP = typename EssentialBC<OpBase>::template Assembly<A>::
361 template BiLinearForm<I>::template OpEssentialLhs<
363
364 auto add_op = [&](auto &bcs) {
366 for (auto &m : bcs) {
367 if (auto bc = m.second->dispBcPtr) {
368 auto &bc_id = m.first;
369 auto regex_str =
370 (boost::format("%s_%s_(.*)") % problem_name % field_name).str();
371 if (std::regex_match(bc_id, std::regex(regex_str))) {
372 MOFEM_TAG_AND_LOG("SELF", Sev::noisy, "OpEssentialLhs") << *bc;
373 pipeline.push_back(
374 new OpSetBc(field_name, false, m.second->getBcMarkersPtr()));
375 pipeline.push_back(new OP(field_name, m.second->getBcEntsPtr()));
376 pipeline.push_back(new OpUnSetBc(field_name));
377 }
378 }
379 }
380 MOFEM_LOG_CHANNEL("SELF");
382 };
383
384 CHKERR add_op(m_field.getInterface<BcManager>()->getBcMapByBlockName());
385
387 }
388};
389
390} // namespace MoFEM
391
392#endif // _ESSENTIAL_DISPLACEMENTCUBITBCDATA_HPP__
static Index< 'o', 3 > o
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
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()
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
#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 double omega
Save field DOFS on vertices/tags.
FTensor::Index< 'm', SPACE_DIM > m
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.
Definition: LogManager.hpp:33
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
constexpr int BASE_DIM
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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:59
constexpr auto field_name
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.
Definition: Essential.hpp:117
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:106
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:200
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
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.
Essential boundary conditions.
Definition: Essential.hpp:125
boost::weak_ptr< FEMethod > fePtr
double vDiag
SmartPetscObj< AO > vAO
SmartPetscObj< Mat > vLhs
MoFEM::Interface & mField
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition: Essential.hpp:47
boost::weak_ptr< FEMethod > fePtr
MoFEM::Interface & mField
double vDiag
SmartPetscObj< Vec > vRhs
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition: Essential.hpp:55
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
Class (Function) to calculate residual side diagonal.
Definition: Essential.hpp:63
typename FormsIntegrators< OpBase >::template Assembly< A >::template BiLinearForm< I >::template OpMass< BASE_DIM, FIELD_DIM > OpMass
Enforce essential constrains on lhs.
Definition: Essential.hpp:95
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:79
Set indices on entities on finite element.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.