v0.13.1
Loading...
Searching...
No Matches
ContactOps.hpp
Go to the documentation of this file.
1
2
3/**
4 * \file ContactOps.hpp
5 * \example ContactOps.hpp
6 */
7
8namespace ContactOps {
9
10
11//! [Common data]
12struct CommonData {
13 boost::shared_ptr<MatrixDouble> mDPtr;
14 boost::shared_ptr<MatrixDouble> mGradPtr;
15 boost::shared_ptr<MatrixDouble> mStrainPtr;
16 boost::shared_ptr<MatrixDouble> mStressPtr;
17
18 boost::shared_ptr<MatrixDouble> contactStressPtr;
19 boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
20 boost::shared_ptr<MatrixDouble> contactTractionPtr;
21 boost::shared_ptr<MatrixDouble> contactDispPtr;
22
23 boost::shared_ptr<MatrixDouble> curlContactStressPtr;
24};
25//! [Common data]
26
31
33 OpConstrainBoundaryRhs(const std::string field_name,
34 boost::shared_ptr<CommonData> common_data_ptr);
35 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
36
37private:
38 boost::shared_ptr<CommonData> commonDataPtr;
39};
40
42 OpConstrainBoundaryLhs_dU(const std::string row_field_name,
43 const std::string col_field_name,
44 boost::shared_ptr<CommonData> common_data_ptr);
45 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
46 EntitiesFieldData::EntData &col_data);
47
48private:
49 boost::shared_ptr<CommonData> commonDataPtr;
50};
51
54 const std::string row_field_name, const std::string col_field_name,
55 boost::shared_ptr<CommonData> common_data_ptr);
56 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
57 EntitiesFieldData::EntData &col_data);
58
59private:
60 boost::shared_ptr<CommonData> commonDataPtr;
61};
62
63template <typename T1, typename T2>
68 t_normal(i) = 0;
69 t_normal(1) = 1.;
70 return t_normal;
71}
72
73template <typename T>
74inline double gap0(FTensor::Tensor1<T, 3> &t_coords,
76 return (-0.5 - t_coords(1)) * t_normal(1);
77}
78
79template <typename T>
80inline double gap(FTensor::Tensor1<T, SPACE_DIM> &t_disp,
82 return t_disp(i) * t_normal(i);
83}
84
85template <typename T>
88 return -t_traction(i) * t_normal(i);
89}
90
91inline double sign(double x) {
92 if (x == 0)
93 return 0;
94 else if (x > 0)
95 return 1;
96 else
97 return -1;
98};
99
100inline double w(const double g, const double t) { return g - cn * t; }
101
102inline double constrian(double &&g0, double &&g, double &&t) {
103 return (w(g - g0, t) + std::abs(w(g - g0, t))) / 2 + g0;
104};
105
106inline double diff_constrains_dtraction(double &&g0, double &&g, double &&t) {
107 return -cn * (1 + sign(w(g - g0, t))) / 2;
108}
109
110inline double diff_constrains_dgap(double &&g0, double &&g, double &&t) {
111 return (1 + sign(w(g - g0, t))) / 2;
112}
113
115 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
117 commonDataPtr(common_data_ptr) {}
118
119MoFEMErrorCode
120OpConstrainBoundaryRhs::iNtegrate(EntitiesFieldData::EntData &data) {
122
123 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
124
125 auto &nf = AssemblyBoundaryEleOp::locF;
126
127 auto t_normal = getFTensor1Normal();
128 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
129
130 auto t_w = getFTensor0IntegrationWeight();
131 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
132 auto t_traction =
133 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
134 auto t_coords = getFTensor1CoordsAtGaussPts();
135
136 size_t nb_base_functions = data.getN().size2() / 3;
137 auto t_base = data.getFTensor1N<3>();
138 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
139
140 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
141
142 const double alpha = t_w * getMeasure();
143
144 auto t_contact_normal = normal(t_coords, t_disp);
146 t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
147
149 t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
150
152
153 t_rhs_constrains(i) =
154 t_contact_normal(i) *
155 constrian(gap0(t_coords, t_contact_normal),
156 gap(t_disp, t_contact_normal),
157 normal_traction(t_traction, t_contact_normal));
158
159 FTensor::Tensor1<double, SPACE_DIM> t_rhs_tangent_disp,
160 t_rhs_tangent_traction;
161 t_rhs_tangent_disp(i) = t_Q(i, j) * t_disp(j);
162 t_rhs_tangent_traction(i) = cn * t_Q(i, j) * t_traction(j);
163
164 size_t bb = 0;
165 for (; bb != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++bb) {
166 const double beta = alpha * (t_base(i) * t_normal(i));
167
168 t_nf(i) -= beta * t_rhs_constrains(i);
169 t_nf(i) -= beta * t_rhs_tangent_disp(i);
170 t_nf(i) += beta * t_rhs_tangent_traction(i);
171
172 ++t_nf;
173 ++t_base;
174 }
175 for (; bb < nb_base_functions; ++bb)
176 ++t_base;
177
178 ++t_disp;
179 ++t_traction;
180 ++t_coords;
181 ++t_w;
182 }
183
185}
186
188 const std::string row_field_name, const std::string col_field_name,
189 boost::shared_ptr<CommonData> common_data_ptr)
190 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
191 DomainEleOp::OPROWCOL),
192 commonDataPtr(common_data_ptr) {
193 sYmm = false;
194}
195
197 EntitiesFieldData::EntData &row_data,
198 EntitiesFieldData::EntData &col_data) {
200
201 const size_t nb_gauss_pts = getGaussPts().size2();
202 auto &locMat = AssemblyBoundaryEleOp::locMat;
203
204 auto t_normal = getFTensor1Normal();
205 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
206
207 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
208 auto t_traction =
209 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
210 auto t_coords = getFTensor1CoordsAtGaussPts();
211
212 auto t_w = getFTensor0IntegrationWeight();
213 auto t_row_base = row_data.getFTensor1N<3>();
214 size_t nb_face_functions = row_data.getN().size2() / 3;
215
216 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
217
218 const double alpha = t_w * getMeasure();
219
220 auto t_contact_normal = normal(t_coords, t_disp);
222 t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
223
225 t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
226
227 auto diff_constrain = diff_constrains_dgap(
228 gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
229 normal_traction(t_traction, t_contact_normal));
230
231 size_t rr = 0;
232 for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
233
234 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
235 locMat, SPACE_DIM * rr);
236
237 const double row_base = t_row_base(i) * t_normal(i);
238
239 auto t_col_base = col_data.getFTensor0N(gg, 0);
240 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
241 ++cc) {
242 const double beta = alpha * row_base * t_col_base;
243
244 t_mat(i, j) -= (beta * diff_constrain) * t_P(i, j);
245 t_mat(i, j) -= beta * t_Q(i, j);
246
247 ++t_col_base;
248 ++t_mat;
249 }
250
251 ++t_row_base;
252 }
253 for (; rr < nb_face_functions; ++rr)
254 ++t_row_base;
255
256 ++t_disp;
257 ++t_traction;
258 ++t_coords;
259 ++t_w;
260 }
261
263}
264
266 const std::string row_field_name, const std::string col_field_name,
267 boost::shared_ptr<CommonData> common_data_ptr)
268 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
269 DomainEleOp::OPROWCOL),
270 commonDataPtr(common_data_ptr) {
271 sYmm = false;
272}
273
275 EntitiesFieldData::EntData &row_data,
276 EntitiesFieldData::EntData &col_data) {
278
279 const size_t nb_gauss_pts = getGaussPts().size2();
280 auto &locMat = AssemblyBoundaryEleOp::locMat;
281
282 auto t_normal = getFTensor1Normal();
283 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
284
285 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
286 auto t_traction =
287 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
288 auto t_coords = getFTensor1CoordsAtGaussPts();
289
290 auto t_w = getFTensor0IntegrationWeight();
291 auto t_row_base = row_data.getFTensor1N<3>();
292 size_t nb_face_functions = row_data.getN().size2() / 3;
293
294 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
295
296 const double alpha = t_w * getMeasure();
297
298 auto t_contact_normal = normal(t_coords, t_disp);
300 t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
301
303 t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
304
305 const double diff_traction = diff_constrains_dtraction(
306 gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
307 normal_traction(t_traction, t_contact_normal));
308
309 size_t rr = 0;
310 for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
311
312 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
313 locMat, SPACE_DIM * rr);
314 const double row_base = t_row_base(i) * t_normal(i);
315
316 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
317 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
318 ++cc) {
319 const double col_base = t_col_base(i) * t_normal(i);
320 const double beta = alpha * row_base * col_base;
321
322 t_mat(i, j) += (beta * diff_traction) * t_P(i, j);
323 t_mat(i, j) += beta * cn * t_Q(i, j);
324
325 ++t_col_base;
326 ++t_mat;
327 }
328
329 ++t_row_base;
330 }
331 for (; rr < nb_face_functions; ++rr)
332 ++t_row_base;
333
334 ++t_disp;
335 ++t_traction;
336 ++t_coords;
337 ++t_w;
338 }
339
341}
342
343}; // namespace ContactOps
constexpr int SPACE_DIM
constexpr double cn
Definition: contact.cpp:124
#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
double normal_traction(FTensor::Tensor1< T, SPACE_DIM > &t_traction, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:86
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:29
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:74
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:28
double w(const double g, const double t)
Definition: ContactOps.hpp:100
double diff_constrains_dgap(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:110
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:65
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:27
double diff_constrains_dtraction(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:106
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:30
double constrian(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:102
double gap(FTensor::Tensor1< T, SPACE_DIM > &t_disp, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:80
double sign(double x)
Definition: ContactOps.hpp:91
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
constexpr double g
boost::shared_ptr< MatrixDouble > curlContactStressPtr
Definition: ContactOps.hpp:23
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ContactOps.hpp:13
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ContactOps.hpp:14
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ContactOps.hpp:15
boost::shared_ptr< MatrixDouble > contactTractionPtr
Definition: ContactOps.hpp:20
boost::shared_ptr< MatrixDouble > contactDispPtr
Definition: ContactOps.hpp:21
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ContactOps.hpp:16
boost::shared_ptr< MatrixDouble > contactStressPtr
Definition: ContactOps.hpp:18
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
Definition: ContactOps.hpp:19
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:265
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:60
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:274
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:49
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:187
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:196
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:114
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:38
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:120