v0.14.0
Loading...
Searching...
No Matches
ElasticSpring.hpp
Go to the documentation of this file.
1/**
2 * @file ElasticSpring.hpp
3 * @author your name (you@domain.com)
4 * @brief Implementation of elastic spring bc
5 * @version 0.13.2
6 * @date 2022-09-18
7 *
8 * @copyright Copyright (c) 2022
9 *
10 */
11
12namespace ElasticExample {
13
14template <CubitBC BC> struct SpringBcType {};
15
16template <CubitBC> struct GetSpringStiffness {
18
19 static MoFEMErrorCode getStiffness(double &normal_stiffness,
20 double &tangent_stiffness,
21 boost::shared_ptr<Range> &ents,
22 MoFEM::Interface &m_field, int ms_id) {
24
25 auto cubit_meshset_ptr =
26 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(ms_id,
27 BLOCKSET);
28 std::vector<double> block_data;
29 CHKERR cubit_meshset_ptr->getAttributes(block_data);
30 if (block_data.size() != 2) {
31 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
32 "Expected that block has two attribute");
33 }
34 normal_stiffness = block_data[0];
35 tangent_stiffness = block_data[1];
36
37 MOFEM_LOG_CHANNEL("WORLD");
38 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "SpringBc")
39 << "Normal stiffness " << normal_stiffness;
40 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "SpringBc")
41 << "Tangent stiffness " << tangent_stiffness;
42
43 ents = boost::make_shared<Range>();
44 CHKERR
45 m_field.get_moab().get_entities_by_handle(cubit_meshset_ptr->meshset,
46 *(ents), true);
47
48 MOFEM_LOG_CHANNEL("SYNC");
49 MOFEM_TAG_AND_LOG("SYNC", Sev::noisy, "SpringBc") << *ents;
50 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::noisy);
51
53 }
54};
55
56} // namespace ElasticExample
57
58template <int FIELD_DIM, AssemblyType A, typename EleOp>
59struct OpFluxRhsImpl<ElasticExample::SpringBcType<BLOCKSET>, 1, FIELD_DIM, A,
60 GAUSS, EleOp> : OpBaseImpl<A, EleOp> {
61
63
64 OpFluxRhsImpl(MoFEM::Interface &m_field, int ms_id, std::string field_name,
65 boost::shared_ptr<MatrixDouble> u_ptr, double scale);
66
67protected:
70 double rhsScale;
71 boost::shared_ptr<MatrixDouble> uPtr;
72 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
73};
74
75template <int FIELD_DIM, AssemblyType A, typename EleOp>
76struct OpFluxLhsImpl<ElasticExample::SpringBcType<BLOCKSET>, 1, FIELD_DIM, A,
77 GAUSS, EleOp> : OpBaseImpl<A, EleOp> {
78
80
81 OpFluxLhsImpl(MoFEM::Interface &m_field, int ms_id,
82 std::string row_field_name, std::string col_field_name);
83
84protected:
87 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
88 EntitiesFieldData::EntData &col_data);
89};
90
91template <int FIELD_DIM, AssemblyType A, typename EleOp>
92OpFluxRhsImpl<ElasticExample::SpringBcType<BLOCKSET>, 1, FIELD_DIM, A, GAUSS,
93 EleOp>::OpFluxRhsImpl(MoFEM::Interface &m_field, int ms_id,
94 std::string field_name,
95 boost::shared_ptr<MatrixDouble> u_ptr,
96 double scale)
97 : OpBase(field_name, field_name, OpBase::OPROW), uPtr(u_ptr), rhsScale(scale) {
99 this->normalStiffness, this->tangentStiffness,
100 this->entsPtr, m_field, ms_id),
101 "Can not read spring stiffness data from blockset");
102}
103
104template <int FIELD_DIM, AssemblyType A, typename EleOp>
105MoFEMErrorCode
107 GAUSS, EleOp>::iNtegrate(EntitiesFieldData::EntData
108 &row_data) {
112 // get element volume
113 const double vol = OpBase::getMeasure();
114 // get integration weights
115 auto t_w = OpBase::getFTensor0IntegrationWeight();
116 // get base function gradient on rows
117 auto t_row_base = row_data.getFTensor0N();
118
119 // get coordinate at integration points
120 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
121
122 // get displacements
123 auto t_u = getFTensor1FromMat<FIELD_DIM>(*uPtr);
124 // loop over integration points
125 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
126 // take into account Jacobian
127 const double alpha = t_w * vol * rhsScale;
128
129 auto l2_norm = t_normal(i) * t_normal(i);
130
131 // normal projection matrix
133 t_P(i, j) = t_normal(i) * t_normal(j) / l2_norm;
134 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
135 // tangential projection matrix
137 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
138 // spring stiffness
140 t_D(i, j) = normalStiffness * t_P(i, j) + tangentStiffness * t_Q(i, j);
141
142 // calculate spring resistance
144 t_reaction(i) = t_D(i, j) * t_u(j);
145
146 // loop over rows base functions
147 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(OpBase::locF);
148 int rr = 0;
149 for (; rr != OpBase::nbRows / FIELD_DIM; ++rr) {
150 t_nf(i) += (alpha * t_row_base) * t_reaction(i);
151 ++t_row_base;
152 ++t_nf;
153 }
154
155 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
156 ++t_row_base;
157
158 ++t_w;
159 ++t_u;
160 ++t_normal;
161 }
163}
164
165template <int FIELD_DIM, AssemblyType A, typename EleOp>
166OpFluxLhsImpl<ElasticExample::SpringBcType<BLOCKSET>, 1, FIELD_DIM, A, GAUSS,
167 EleOp>::OpFluxLhsImpl(MoFEM::Interface &m_field, int ms_id,
168 std::string row_field_name,
169 std::string col_field_name)
170 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL) {
172 this->normalStiffness, this->tangentStiffness,
173 this->entsPtr, m_field, ms_id),
174 "Can not read spring stiffness data from blockset");
175}
176
177template <int FIELD_DIM, AssemblyType A, typename EleOp>
178MoFEMErrorCode
180 GAUSS,
181 EleOp>::iNtegrate(EntitiesFieldData::EntData &row_data,
182 EntitiesFieldData::EntData &col_data) {
183
186
188 // get element volume
189 const double vol = OpBase::getMeasure();
190 // get integration weights
191 auto t_w = OpBase::getFTensor0IntegrationWeight();
192 // get base function gradient on rows
193 auto t_row_base = row_data.getFTensor0N();
194
195 // get coordinate at integration points
196 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
197
198 // loop over integration points
199 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
200 // take into account Jacobian
201 const double alpha = t_w * vol;
202
203 auto l2_norm = t_normal(i) * t_normal(i);
204
205 // normal projection matrix
207 t_P(i, j) = t_normal(i) * t_normal(j) / l2_norm;
208 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
209 // tangential projection matrix
211 t_Q(i, j) = t_kd(i, j) - t_P(i, j);
212 // spring stiffness
214 t_D(i, j) = normalStiffness * t_P(i, j) + tangentStiffness * t_Q(i, j);
215
216 // loop over rows base functions
217 int rr = 0;
218 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
219 // get column base functions gradient at gauss point gg
220 auto t_col_base = col_data.getFTensor0N(gg, 0);
221 // get sub matrix for the row
222 auto t_m = OpBase::template getLocMat<FIELD_DIM>(FIELD_DIM * rr);
223 // loop over columns
224 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
225 // calculate element of local matrix
226 t_m(i, j) += t_D(i, j) * (alpha * t_row_base * t_col_base);
227 ++t_col_base;
228 ++t_m;
229 }
230 ++t_row_base;
231 }
232 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
233 ++t_row_base;
234
235 ++t_normal;
236 ++t_w; // move to another integration weight
237 }
239}
240
241template <CubitBC BCTYPE, int BASE_DIM, int FIELD_DIM, AssemblyType A,
242 IntegrationType I, typename OpBase>
243struct AddFluxToRhsPipelineImpl<
244
245 OpFluxRhsImpl<ElasticExample::SpringBcType<BCTYPE>, BASE_DIM, FIELD_DIM, A,
246 I, OpBase>,
247 A, I, OpBase
248
249 > {
250
252
253 static MoFEMErrorCode add(
254
255 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
256 MoFEM::Interface &m_field, std::string field_name,
257 boost::shared_ptr<MatrixDouble> u_ptr, double scale,
258 std::string block_name, Sev sev
259
260 ) {
262
263 using OP = OpFluxRhsImpl<ElasticExample::SpringBcType<BLOCKSET>, BASE_DIM,
265
266 auto add_op = [&](auto &&meshset_vec_ptr) {
267 for (auto m : meshset_vec_ptr) {
268 MOFEM_TAG_AND_LOG("WORLD", sev, "OpSpringRhs") << "Add " << *m;
269 pipeline.push_back(
270 new OP(m_field, m->getMeshsetId(), field_name, u_ptr, scale));
271 }
272 MOFEM_LOG_CHANNEL("WORLD");
273 };
274
275 switch (BCTYPE) {
276 case BLOCKSET:
277 add_op(
278
279 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
280 std::regex(
281
282 (boost::format("%s(.*)") % block_name).str()
283
284 ))
285
286 );
287
288 break;
289 default:
290 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
291 "Handling of bc type not implemented");
292 break;
293 }
295 }
296};
297
298template <CubitBC BCTYPE, int BASE_DIM, int FIELD_DIM, AssemblyType A,
299 IntegrationType I, typename OpBase>
300struct AddFluxToLhsPipelineImpl<
301
302 OpFluxLhsImpl<ElasticExample::SpringBcType<BCTYPE>, BASE_DIM, FIELD_DIM, A,
303 I, OpBase>,
304 A, I, OpBase
305
306 > {
307
309
310 static MoFEMErrorCode add(
311
312 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
313 MoFEM::Interface &m_field, std::string row_field_name,
314 std::string col_field_name, std::string block_name, Sev sev
315
316 ) {
318
319 using OP = OpFluxLhsImpl<ElasticExample::SpringBcType<BLOCKSET>, BASE_DIM,
321
322 auto add_op = [&](auto &&meshset_vec_ptr) {
323 for (auto m : meshset_vec_ptr) {
324 MOFEM_TAG_AND_LOG("WORLD", sev, "OpSpringLhs") << "Add " << *m;
325 pipeline.push_back(
326 new OP(m_field, m->getMeshsetId(), row_field_name, col_field_name));
327 }
328 MOFEM_LOG_CHANNEL("WORLD");
329 };
330
331 switch (BCTYPE) {
332 case BLOCKSET:
333 add_op(
334
335 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
336 std::regex(
337
338 (boost::format("%s(.*)") % block_name).str()
339
340 ))
341
342 );
343
344 break;
345 default:
346 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
347 "Handling of bc type not implemented");
348 break;
349 }
351 }
352};
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition: LogManager.hpp:352
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
constexpr int FIELD_DIM
Kronecker Delta class.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
CubitBC
Types of sets and boundary conditions.
Definition: definitions.h:144
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
constexpr int BASE_DIM
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
constexpr IntegrationType I
constexpr AssemblyType A
double scale
Definition: plastic.cpp:166
constexpr auto field_name
static MoFEMErrorCode add(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, MoFEM::Interface &m_field, std::string row_field_name, std::string col_field_name, std::string block_name, Sev sev)
static MoFEMErrorCode add(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, MoFEM::Interface &m_field, std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, double scale, std::string block_name, Sev sev)
static MoFEMErrorCode getStiffness(double &normal_stiffness, double &tangent_stiffness, boost::shared_ptr< Range > &ents, MoFEM::Interface &m_field, int ms_id)
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)