v0.10.0
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 namespace ContactOps {
16 
17 //! [Common data]
18 struct CommonData {
19  boost::shared_ptr<MatrixDouble> mDPtr;
20  boost::shared_ptr<MatrixDouble> mGradPtr;
21  boost::shared_ptr<MatrixDouble> mStrainPtr;
22  boost::shared_ptr<MatrixDouble> mStressPtr;
23 
24  boost::shared_ptr<MatrixDouble> contactStressPtr;
25  boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
26  boost::shared_ptr<MatrixDouble> contactTractionPtr;
27  boost::shared_ptr<MatrixDouble> contactDispPtr;
28 
29  boost::shared_ptr<MatrixDouble> curlContactStressPtr;
30 };
31 //! [Common data]
32 
37 
39  OpConstrainBoundaryRhs(const std::string field_name,
40  boost::shared_ptr<CommonData> common_data_ptr);
41  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
42 
43 private:
44  boost::shared_ptr<CommonData> commonDataPtr;
45 };
46 
48  OpConstrainBoundaryLhs_dU(const std::string row_field_name,
49  const std::string col_field_name,
50  boost::shared_ptr<CommonData> common_data_ptr);
51  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
52  EntityType col_type, EntData &row_data,
53  EntData &col_data);
54 
55 private:
56  boost::shared_ptr<CommonData> commonDataPtr;
58 };
59 
62  const std::string row_field_name, const std::string col_field_name,
63  boost::shared_ptr<CommonData> common_data_ptr);
64  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
65  EntityType col_type, EntData &row_data,
66  EntData &col_data);
67 
68 private:
69  boost::shared_ptr<CommonData> commonDataPtr;
71 };
72 
73 template <typename T1, typename T2>
78  t_normal(i) = 0;
79  t_normal(1) = 1.;
80  return t_normal;
81 }
82 
83 template <typename T>
84 inline double gap0(FTensor::Tensor1<T, 3> &t_coords,
86  return (-0.5 - t_coords(1)) * t_normal(1);
87 }
88 
89 template <typename T>
90 inline double gap(FTensor::Tensor1<T, SPACE_DIM> &t_disp,
92  return t_disp(i) * t_normal(i);
93 }
94 
95 template <typename T>
98  return -t_traction(i) * t_normal(i);
99 }
100 
101 inline double sign(double x) {
102  if (x == 0)
103  return 0;
104  else if (x > 0)
105  return 1;
106  else
107  return -1;
108 };
109 
110 inline double w(const double g, const double t) { return g - cn * t; }
111 
112 inline double constrian(double &&g0, double &&g, double &&t) {
113  return (w(g - g0, t) + std::abs(w(g - g0, t))) / 2 + g0;
114 };
115 
116 inline double diff_constrains_dtraction(double &&g0, double &&g, double &&t) {
117  return -cn * (1 + sign(w(g - g0, t))) / 2;
118 }
119 
120 inline double diff_constrains_dgap(double &&g0, double &&g, double &&t) {
121  return (1 + sign(w(g - g0, t))) / 2;
122 }
123 
124 
126  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
127  : BoundaryEleOp(field_name, DomainEleOp::OPROW),
128  commonDataPtr(common_data_ptr) {}
129 
131  EntData &data) {
133 
134  const size_t nb_gauss_pts = getGaussPts().size2();
135  const size_t nb_dofs = data.getIndices().size();
136 
137  if (nb_dofs) {
138 
139  std::array<double, MAX_DOFS_ON_ENTITY> nf;
140  std::fill(&nf[0], &nf[nb_dofs], 0);
141 
142  auto t_normal = getFTensor1Normal();
143  t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
144 
145  auto t_w = getFTensor0IntegrationWeight();
146  auto t_disp =
147  getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
148  auto t_traction =
149  getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
150  auto t_coords = getFTensor1CoordsAtGaussPts();
151 
152  size_t nb_base_functions = data.getN().size2() / 3;
153  auto t_base = data.getFTensor1N<3>();
154  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
155 
156  auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data());
157 
158  const double alpha = t_w * getMeasure();
159 
160  auto t_contact_normal = normal(t_coords, t_disp);
162  t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
163 
165  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
166 
167  FTensor::Tensor1<double, SPACE_DIM> t_rhs_constrains;
168 
169  t_rhs_constrains(i) =
170  t_contact_normal(i) *
171  constrian(gap0(t_coords, t_contact_normal),
172  gap(t_disp, t_contact_normal),
173  normal_traction(t_traction, t_contact_normal));
174 
175  FTensor::Tensor1<double, SPACE_DIM> t_rhs_tangent_disp,
176  t_rhs_tangent_traction;
177  t_rhs_tangent_disp(i) = t_Q(i, j) * t_disp(j);
178  t_rhs_tangent_traction(i) = cn * t_Q(i, j) * t_traction(j);
179 
180  size_t bb = 0;
181  for (; bb != nb_dofs / SPACE_DIM; ++bb) {
182  const double beta = alpha * (t_base(i) * t_normal(i));
183 
184  t_nf(i) -= beta * t_rhs_constrains(i);
185  t_nf(i) -= beta * t_rhs_tangent_disp(i);
186  t_nf(i) += beta * t_rhs_tangent_traction(i);
187 
188  ++t_nf;
189  ++t_base;
190  }
191  for (; bb < nb_base_functions; ++bb)
192  ++t_base;
193 
194  ++t_disp;
195  ++t_traction;
196  ++t_coords;
197  ++t_w;
198  }
199 
200  CHKERR VecSetValues(getSNESf(), data, nf.data(), ADD_VALUES);
201  }
202 
204 }
205 
207  const std::string row_field_name, const std::string col_field_name,
208  boost::shared_ptr<CommonData> common_data_ptr)
209  : BoundaryEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
210  commonDataPtr(common_data_ptr) {
211  sYmm = false;
212 }
214  EntityType row_type,
215  EntityType col_type,
216  EntData &row_data,
217  EntData &col_data) {
219 
220  const size_t nb_gauss_pts = getGaussPts().size2();
221  const size_t row_nb_dofs = row_data.getIndices().size();
222  const size_t col_nb_dofs = col_data.getIndices().size();
223 
224  if (row_nb_dofs && col_nb_dofs) {
225 
226  locMat.resize(row_nb_dofs, col_nb_dofs, false);
227  locMat.clear();
228 
229  auto t_normal = getFTensor1Normal();
230  t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
231 
232  auto t_disp =
233  getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
234  auto t_traction =
235  getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
236  auto t_coords = getFTensor1CoordsAtGaussPts();
237 
238  auto t_w = getFTensor0IntegrationWeight();
239  auto t_row_base = row_data.getFTensor1N<3>();
240  size_t nb_face_functions = row_data.getN().size2() / 3;
241 
242  locMat.resize(row_nb_dofs, col_nb_dofs, false);
243  locMat.clear();
244 
245  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
246 
247  const double alpha = t_w * getMeasure();
248 
249  auto t_contact_normal = normal(t_coords, t_disp);
251  t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
252 
254  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
255 
256  auto diff_constrain = diff_constrains_dgap(
257  gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
258  normal_traction(t_traction, t_contact_normal));
259 
260  size_t rr = 0;
261  for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
262 
263  auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
264  locMat, SPACE_DIM * rr);
265 
266  const double row_base = t_row_base(i) * t_normal(i);
267 
268  auto t_col_base = col_data.getFTensor0N(gg, 0);
269  for (size_t cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
270  const double beta = alpha * row_base * t_col_base;
271 
272  t_mat(i, j) -= (beta * diff_constrain) * t_P(i, j);
273  t_mat(i, j) -= beta * t_Q(i, j);
274 
275  ++t_col_base;
276  ++t_mat;
277  }
278 
279  ++t_row_base;
280  }
281  for (; rr < nb_face_functions; ++rr)
282  ++t_row_base;
283 
284  ++t_disp;
285  ++t_traction;
286  ++t_coords;
287  ++t_w;
288  }
289 
290  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
291  ADD_VALUES);
292  }
293 
295 }
296 
298  const std::string row_field_name, const std::string col_field_name,
299  boost::shared_ptr<CommonData> common_data_ptr)
300  : BoundaryEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
301  commonDataPtr(common_data_ptr) {
302  sYmm = false;
303 }
305  int row_side, int col_side, EntityType row_type, EntityType col_type,
306  EntData &row_data, EntData &col_data) {
308 
309  const size_t nb_gauss_pts = getGaussPts().size2();
310  const size_t row_nb_dofs = row_data.getIndices().size();
311  const size_t col_nb_dofs = col_data.getIndices().size();
312 
313  if (row_nb_dofs && col_nb_dofs) {
314 
315  locMat.resize(row_nb_dofs, col_nb_dofs, false);
316  locMat.clear();
317 
318  auto t_normal = getFTensor1Normal();
319  t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
320 
321  auto t_disp =
322  getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
323  auto t_traction =
324  getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
325  auto t_coords = getFTensor1CoordsAtGaussPts();
326 
327  auto t_w = getFTensor0IntegrationWeight();
328  auto t_row_base = row_data.getFTensor1N<3>();
329  size_t nb_face_functions = row_data.getN().size2() / 3;
330 
331  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
332 
333  const double alpha = t_w * getMeasure();
334 
335  auto t_contact_normal = normal(t_coords, t_disp);
337  t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
338 
340  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
341 
342  const double diff_traction = diff_constrains_dtraction(
343  gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
344  normal_traction(t_traction, t_contact_normal));
345 
346  size_t rr = 0;
347  for (; rr != row_nb_dofs / SPACE_DIM; ++rr) {
348 
349  auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
350  locMat, SPACE_DIM * rr);
351  const double row_base = t_row_base(i) * t_normal(i);
352 
353  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
354  for (size_t cc = 0; cc != col_nb_dofs / SPACE_DIM; ++cc) {
355  const double col_base = t_col_base(i) * t_normal(i);
356  const double beta = alpha * row_base * col_base;
357 
358  t_mat(i, j) += (beta * diff_traction) * t_P(i, j);
359  t_mat(i, j) += beta * cn * t_Q(i, j);
360 
361  ++t_col_base;
362  ++t_mat;
363  }
364 
365  ++t_row_base;
366  }
367  for (; rr < nb_face_functions; ++rr)
368  ++t_row_base;
369 
370  ++t_disp;
371  ++t_traction;
372  ++t_coords;
373  ++t_w;
374  }
375 
376  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
377  ADD_VALUES);
378  }
379 
381 }
382 
383 }; // namespace OpContactTools
ContactOps
Definition: ContactOps.hpp:15
ContactOps::k
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:35
cn
constexpr double cn
Definition: plastic.cpp:104
FTensor::Tensor2
Definition: Tensor2_value.hpp:17
ContactOps::normal_traction
double normal_traction(FTensor::Tensor1< T, SPACE_DIM > &t_traction, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:96
ContactOps::CommonData::contactTractionPtr
boost::shared_ptr< MatrixDouble > contactTractionPtr
Definition: ContactOps.hpp:26
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1020
ContactOps::OpConstrainBoundaryLhs_dU::OpConstrainBoundaryLhs_dU
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:206
MoFEM::ForcesAndSourcesCore::UserDataOperator::getSNESf
Vec getSNESf() const
Definition: ForcesAndSourcesCore.hpp:960
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM::DataForcesAndSourcesCore::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: DataStructures.hpp:1178
MoFEM::Types::MatrixDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
ContactOps::OpConstrainBoundaryLhs_dTraction::locMat
MatrixDouble locMat
Definition: ContactOps.hpp:70
FTensor::Tensor1
Definition: Tensor1_value.hpp:9
FTensor::kronecker_delta
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
Definition: Kronecker_Delta.hpp:81
ContactOps::w
double w(const double g, const double t)
Definition: ContactOps.hpp:110
ContactOps::CommonData::mStrainPtr
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ContactOps.hpp:21
ContactOps::CommonData::contactStressDivergencePtr
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
Definition: ContactOps.hpp:25
MoFEM::FaceElementForcesAndSourcesCoreBase::UserDataOperator::getMeasure
double getMeasure()
get measure of element
Definition: FaceElementForcesAndSourcesCore.hpp:419
ContactOps::OpConstrainBoundaryLhs_dU
Definition: ContactOps.hpp:47
MoFEM::FaceElementForcesAndSourcesCoreBase::UserDataOperator::getFTensor1Normal
auto getFTensor1Normal()
get normal as tensor
Definition: FaceElementForcesAndSourcesCore.hpp:439
ContactOps::diff_constrains_dgap
double diff_constrains_dgap(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:120
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition: DataStructures.cpp:587
ContactOps::OpConstrainBoundaryLhs_dU::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:56
ElectroPhysiology::alpha
const double alpha
Definition: ElecPhysOperators.hpp:28
SPACE_DIM
constexpr int SPACE_DIM
Definition: plastic.cpp:50
ContactOps::OpConstrainBoundaryLhs_dU::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: ContactOps.hpp:213
convert.type
type
Definition: convert.py:66
ContactOps::constrian
double constrian(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:112
MoFEM::FaceElementForcesAndSourcesCoreBase::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:44
Poisson2DTraditionalDDOperators::MatSetValues
constexpr auto MatSetValues
Definition: poisson_dd_H.hpp:20
ContactOps::sign
double sign(double x)
Definition: ContactOps.hpp:101
FTensor::Index< 'i', SPACE_DIM >
ContactOps::OpConstrainBoundaryRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: ContactOps.hpp:130
ContactOps::OpConstrainBoundaryLhs_dU::locMat
MatrixDouble locMat
Definition: ContactOps.hpp:57
ContactOps::CommonData::mGradPtr
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ContactOps.hpp:20
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: DataStructures.hpp:1459
ContactOps::diff_constrains_dtraction
double diff_constrains_dtraction(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:116
ContactOps::CommonData::contactDispPtr
boost::shared_ptr< MatrixDouble > contactDispPtr
Definition: ContactOps.hpp:27
ContactOps::l
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:36
ContactOps::CommonData::mStressPtr
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ContactOps.hpp:22
ContactOps::gap
double gap(FTensor::Tensor1< T, SPACE_DIM > &t_disp, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:90
ContactOps::CommonData
[Common data]
Definition: ContactOps.hpp:18
ContactOps::OpConstrainBoundaryLhs_dTraction
Definition: ContactOps.hpp:60
ContactOps::CommonData::contactStressPtr
boost::shared_ptr< MatrixDouble > contactStressPtr
Definition: ContactOps.hpp:24
ContactOps::OpConstrainBoundaryRhs::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:44
ContactOps::gap0
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, SPACE_DIM > &t_normal)
Definition: ContactOps.hpp:84
ContactOps::OpConstrainBoundaryRhs
Definition: ContactOps.hpp:38
MoFEM::ForcesAndSourcesCore::UserDataOperator::getSNESB
Mat getSNESB() const
Definition: ForcesAndSourcesCore.hpp:972
ContactOps::OpConstrainBoundaryLhs_dTraction::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: ContactOps.hpp:304
ContactOps::CommonData::curlContactStressPtr
boost::shared_ptr< MatrixDouble > curlContactStressPtr
Definition: ContactOps.hpp:29
ContactOps::CommonData::mDPtr
boost::shared_ptr< MatrixDouble > mDPtr
Definition: ContactOps.hpp:19
ContactOps::normal
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:75
MoFEM::DataForcesAndSourcesCore::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: DataStructures.hpp:1288
Poisson2DTraditionalDDOperators::VecSetValues
constexpr auto VecSetValues
Definition: poisson_dd_H.hpp:19
MoFEM::FaceElementForcesAndSourcesCoreBase::UserDataOperator::getFTensor1CoordsAtGaussPts
auto getFTensor1CoordsAtGaussPts()
get coordinates at Gauss pts.
Definition: FaceElementForcesAndSourcesCore.hpp:484
ContactOps::OpConstrainBoundaryLhs_dTraction::OpConstrainBoundaryLhs_dTraction
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:297
MoFEM::DataForcesAndSourcesCore::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: DataStructures.hpp:60
ContactOps::OpConstrainBoundaryLhs_dTraction::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:69
MoFEM::VolumeElementForcesAndSourcesCoreBase::UserDataOperator
default operator for TET element
Definition: VolumeElementForcesAndSourcesCore.hpp:45
ContactOps::i
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:33
ContactOps::OpConstrainBoundaryRhs::OpConstrainBoundaryRhs
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:125
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:72
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ContactOps::j
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:34
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1016