v0.13.1
ContactOps.hpp
Go to the documentation of this file.
1/* This file is part of MoFEM.
2 * MoFEM is free software: you can redistribute it and/or modify it under
3 * the terms of the GNU Lesser General Public License as published by the
4 * Free Software Foundation, either version 3 of the License, or (at your
5 * option) any later version.
6 *
7 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10 * License for more details.
11 *
12 * You should have received a copy of the GNU Lesser General Public
13 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14
15/**
16 * \file ContactOps.hpp
17 * \example ContactOps.hpp
18 */
19
20namespace ContactOps {
21
22
23//! [Common data]
24struct CommonData {
25 boost::shared_ptr<MatrixDouble> mDPtr;
26 boost::shared_ptr<MatrixDouble> mGradPtr;
27 boost::shared_ptr<MatrixDouble> mStrainPtr;
28 boost::shared_ptr<MatrixDouble> mStressPtr;
29
30 boost::shared_ptr<MatrixDouble> contactStressPtr;
31 boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
32 boost::shared_ptr<MatrixDouble> contactTractionPtr;
33 boost::shared_ptr<MatrixDouble> contactDispPtr;
34
35 boost::shared_ptr<MatrixDouble> curlContactStressPtr;
36};
37//! [Common data]
38
43
45 OpConstrainBoundaryRhs(const std::string field_name,
46 boost::shared_ptr<CommonData> common_data_ptr);
48
49private:
50 boost::shared_ptr<CommonData> commonDataPtr;
51};
52
54 OpConstrainBoundaryLhs_dU(const std::string row_field_name,
55 const std::string col_field_name,
56 boost::shared_ptr<CommonData> common_data_ptr);
59
60private:
61 boost::shared_ptr<CommonData> commonDataPtr;
62};
63
66 const std::string row_field_name, const std::string col_field_name,
67 boost::shared_ptr<CommonData> common_data_ptr);
70
71private:
72 boost::shared_ptr<CommonData> commonDataPtr;
73};
74
75template <typename T1, typename T2>
80 t_normal(i) = 0;
81 t_normal(1) = 1.;
82 return t_normal;
83}
84
85template <typename T>
86inline double gap0(FTensor::Tensor1<T, 3> &t_coords,
88 return (-0.5 - t_coords(1)) * t_normal(1);
89}
90
91template <typename T>
92inline double gap(FTensor::Tensor1<T, SPACE_DIM> &t_disp,
94 return t_disp(i) * t_normal(i);
95}
96
97template <typename T>
100 return -t_traction(i) * t_normal(i);
101}
102
103inline double sign(double x) {
104 if (x == 0)
105 return 0;
106 else if (x > 0)
107 return 1;
108 else
109 return -1;
110};
111
112inline double w(const double g, const double t) { return g - cn * t; }
113
114inline double constrian(double &&g0, double &&g, double &&t) {
115 return (w(g - g0, t) + std::abs(w(g - g0, t))) / 2 + g0;
116};
117
118inline double diff_constrains_dtraction(double &&g0, double &&g, double &&t) {
119 return -cn * (1 + sign(w(g - g0, t))) / 2;
120}
121
122inline double diff_constrains_dgap(double &&g0, double &&g, double &&t) {
123 return (1 + sign(w(g - g0, t))) / 2;
124}
125
127 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
129 commonDataPtr(common_data_ptr) {}
130
134
135 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
136
137 auto &nf = AssemblyBoundaryEleOp::locF;
138
139 auto t_normal = getFTensor1Normal();
140 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
141
142 auto t_w = getFTensor0IntegrationWeight();
143 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
144 auto t_traction =
145 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
146 auto t_coords = getFTensor1CoordsAtGaussPts();
147
148 size_t nb_base_functions = data.getN().size2() / 3;
149 auto t_base = data.getFTensor1N<3>();
150 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
151
152 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
153
154 const double alpha = t_w * getMeasure();
155
156 auto t_contact_normal = normal(t_coords, t_disp);
158 t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
159
161 t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
162
164
165 t_rhs_constrains(i) =
166 t_contact_normal(i) *
167 constrian(gap0(t_coords, t_contact_normal),
168 gap(t_disp, t_contact_normal),
169 normal_traction(t_traction, t_contact_normal));
170
171 FTensor::Tensor1<double, SPACE_DIM> t_rhs_tangent_disp,
172 t_rhs_tangent_traction;
173 t_rhs_tangent_disp(i) = t_Q(i, j) * t_disp(j);
174 t_rhs_tangent_traction(i) = cn * t_Q(i, j) * t_traction(j);
175
176 size_t bb = 0;
177 for (; bb != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++bb) {
178 const double beta = alpha * (t_base(i) * t_normal(i));
179
180 t_nf(i) -= beta * t_rhs_constrains(i);
181 t_nf(i) -= beta * t_rhs_tangent_disp(i);
182 t_nf(i) += beta * t_rhs_tangent_traction(i);
183
184 ++t_nf;
185 ++t_base;
186 }
187 for (; bb < nb_base_functions; ++bb)
188 ++t_base;
189
190 ++t_disp;
191 ++t_traction;
192 ++t_coords;
193 ++t_w;
194 }
195
197}
198
200 const std::string row_field_name, const std::string col_field_name,
201 boost::shared_ptr<CommonData> common_data_ptr)
202 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
203 DomainEleOp::OPROWCOL),
204 commonDataPtr(common_data_ptr) {
205 sYmm = false;
206}
207
210 EntitiesFieldData::EntData &col_data) {
212
213 const size_t nb_gauss_pts = getGaussPts().size2();
214 auto &locMat = AssemblyBoundaryEleOp::locMat;
215
216 auto t_normal = getFTensor1Normal();
217 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
218
219 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
220 auto t_traction =
221 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
222 auto t_coords = getFTensor1CoordsAtGaussPts();
223
224 auto t_w = getFTensor0IntegrationWeight();
225 auto t_row_base = row_data.getFTensor1N<3>();
226 size_t nb_face_functions = row_data.getN().size2() / 3;
227
228 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
229
230 const double alpha = t_w * getMeasure();
231
232 auto t_contact_normal = normal(t_coords, t_disp);
234 t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
235
237 t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
238
239 auto diff_constrain = diff_constrains_dgap(
240 gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
241 normal_traction(t_traction, t_contact_normal));
242
243 size_t rr = 0;
244 for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
245
246 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
247 locMat, SPACE_DIM * rr);
248
249 const double row_base = t_row_base(i) * t_normal(i);
250
251 auto t_col_base = col_data.getFTensor0N(gg, 0);
252 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
253 ++cc) {
254 const double beta = alpha * row_base * t_col_base;
255
256 t_mat(i, j) -= (beta * diff_constrain) * t_P(i, j);
257 t_mat(i, j) -= beta * t_Q(i, j);
258
259 ++t_col_base;
260 ++t_mat;
261 }
262
263 ++t_row_base;
264 }
265 for (; rr < nb_face_functions; ++rr)
266 ++t_row_base;
267
268 ++t_disp;
269 ++t_traction;
270 ++t_coords;
271 ++t_w;
272 }
273
275}
276
278 const std::string row_field_name, const std::string col_field_name,
279 boost::shared_ptr<CommonData> common_data_ptr)
280 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
281 DomainEleOp::OPROWCOL),
282 commonDataPtr(common_data_ptr) {
283 sYmm = false;
284}
285
288 EntitiesFieldData::EntData &col_data) {
290
291 const size_t nb_gauss_pts = getGaussPts().size2();
292 auto &locMat = AssemblyBoundaryEleOp::locMat;
293
294 auto t_normal = getFTensor1Normal();
295 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
296
297 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
298 auto t_traction =
299 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
300 auto t_coords = getFTensor1CoordsAtGaussPts();
301
302 auto t_w = getFTensor0IntegrationWeight();
303 auto t_row_base = row_data.getFTensor1N<3>();
304 size_t nb_face_functions = row_data.getN().size2() / 3;
305
306 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
307
308 const double alpha = t_w * getMeasure();
309
310 auto t_contact_normal = normal(t_coords, t_disp);
312 t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
313
315 t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
316
317 const double diff_traction = diff_constrains_dtraction(
318 gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
319 normal_traction(t_traction, t_contact_normal));
320
321 size_t rr = 0;
322 for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
323
324 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
325 locMat, SPACE_DIM * rr);
326 const double row_base = t_row_base(i) * t_normal(i);
327
328 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
329 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
330 ++cc) {
331 const double col_base = t_col_base(i) * t_normal(i);
332 const double beta = alpha * row_base * col_base;
333
334 t_mat(i, j) += (beta * diff_traction) * t_P(i, j);
335 t_mat(i, j) += beta * cn * t_Q(i, j);
336
337 ++t_col_base;
338 ++t_mat;
339 }
340
341 ++t_row_base;
342 }
343 for (; rr < nb_face_functions; ++rr)
344 ++t_row_base;
345
346 ++t_disp;
347 ++t_traction;
348 ++t_coords;
349 ++t_w;
350 }
351
353}
354
355}; // namespace ContactOps
EntitiesFieldData::EntData EntData
constexpr int SPACE_DIM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
constexpr double beta
constexpr double alpha
double normal_traction(FTensor::Tensor1< T, SPACE_DIM > &t_traction, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:98
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:41
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:86
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:40
double w(const double g, const double t)
Definition: ContactOps.hpp:112
double diff_constrains_dgap(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:122
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:39
double diff_constrains_dtraction(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:118
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:42
double constrian(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:114
double gap(FTensor::Tensor1< T, SPACE_DIM > &t_disp, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:92
double sign(double x)
Definition: ContactOps.hpp:103
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
double cn
Definition: plastic.cpp:111
constexpr double t
plate stiffness
Definition: plate.cpp:76
constexpr auto field_name
constexpr double g
boost::shared_ptr< MatrixDouble > curlContactStressPtr
Definition: ContactOps.hpp:35
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ContactOps.hpp:25
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ContactOps.hpp:26
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ContactOps.hpp:27
boost::shared_ptr< MatrixDouble > contactTractionPtr
Definition: ContactOps.hpp:32
boost::shared_ptr< MatrixDouble > contactDispPtr
Definition: ContactOps.hpp:33
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ContactOps.hpp:28
boost::shared_ptr< MatrixDouble > contactStressPtr
Definition: ContactOps.hpp:30
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
Definition: ContactOps.hpp:31
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:277
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:72
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:286
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:61
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:199
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:208
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:126
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:50
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:132