v0.16.0
Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | List of all members
OpCalculateRotationAndSpatialGradient Struct Reference

#include "users_modules/eshelbian_plasticity/src/EshelbianOperators.hpp"

Inheritance diagram for OpCalculateRotationAndSpatialGradient:
[legend]
Collaboration diagram for OpCalculateRotationAndSpatialGradient:
[legend]

Public Member Functions

 OpCalculateRotationAndSpatialGradient (boost::shared_ptr< DataAtIntegrationPts > data_ptr)
 
MoFEMErrorCode doWork (int side, EntityType type, EntData &data)
 Operator for linear form, usually to calculate values on right hand side.
 
- Public Member Functions inherited from MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
int getNumNodes ()
 get element number of nodes
 
const EntityHandle * getConn ()
 get element connectivity
 
double getVolume () const
 element volume (linear geometry)
 
double & getVolume ()
 element volume (linear geometry)
 
FTensor::Tensor2< double *, 3, 3 > & getJac ()
 get element Jacobian
 
FTensor::Tensor2< double *, 3, 3 > & getInvJac ()
 get element inverse Jacobian
 
VectorDouble & getCoords ()
 nodal coordinates
 
VolumeElementForcesAndSourcesCore * getVolumeFE () const
 return pointer to Generic Volume Finite Element object
 
- Public Member Functions inherited from MoFEM::ForcesAndSourcesCore::UserDataOperator
 UserDataOperator (const FieldSpace space, const char type=OPSPACE, const bool symm=true)
 Constructor for operators working on finite element spaces.
 
 UserDataOperator (const std::string field_name, const char type, const bool symm=true)
 Constructor for operators working on a single field.
 
 UserDataOperator (const std::string row_field_name, const std::string col_field_name, const char type, const bool symm=true)
 Constructor for operators working on two fields (bilinear forms)
 
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr () const
 Return raw pointer to NumeredEntFiniteElement.
 
EntityHandle getFEEntityHandle () const
 Return finite element entity handle.
 
int getFEDim () const
 Get dimension of finite element.
 
EntityType getFEType () const
 Get dimension of finite element.
 
boost::weak_ptr< SideNumber > getSideNumberPtr (const int side_number, const EntityType type)
 Get the side number pointer.
 
EntityHandle getSideEntity (const int side_number, const EntityType type)
 Get the side entity.
 
int getNumberOfNodesOnElement () const
 Get the number of nodes on finite element.
 
MoFEMErrorCode getProblemRowIndices (const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
 Get row indices.
 
MoFEMErrorCode getProblemColIndices (const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
 Get col indices.
 
const FEMethod * getFEMethod () const
 Return raw pointer to Finite Element Method object.
 
int getOpType () const
 Get operator types.
 
void setOpType (const OpType type)
 Set operator type.
 
void addOpType (const OpType type)
 Add operator type.
 
int getNinTheLoop () const
 get number of finite element in the loop
 
int getLoopSize () const
 get size of elements in the loop
 
std::string getFEName () const
 Get name of the element.
 
ForcesAndSourcesCore * getPtrFE () const
 
ForcesAndSourcesCore * getSidePtrFE () const
 
ForcesAndSourcesCore * getRefinePtrFE () const
 
const PetscData::Switches & getDataCtx () const
 
KspMethod::KSPContext getKSPCtx () const
 
SnesMethod::SNESContext getSNESCtx () const
 
TSMethod::TSContext getTSCtx () const
 
Vec getKSPf () const
 
Mat getKSPA () const
 
Mat getKSPB () const
 
Vec getSNESf () const
 
Vec getSNESx () const
 
Mat getSNESA () const
 
Mat getSNESB () const
 
Vec getTSu () const
 
Vec getTSu_t () const
 
Vec getTSu_tt () const
 
Vec getTSf () const
 
Mat getTSA () const
 
Mat getTSB () const
 
int getTSstep () const
 
double getTStime () const
 
double getTStimeStep () const
 
double getTSa () const
 
double getTSaa () const
 
MatrixDouble & getGaussPts ()
 matrix of integration (Gauss) points for Volume Element
 
auto getFTensor0IntegrationWeight ()
 Get integration weights.
 
MatrixDouble & getCoordsAtGaussPts ()
 Gauss points and weight, matrix (nb. of points x 3)
 
auto getFTensor1CoordsAtGaussPts ()
 Get coordinates at integration points assuming linear geometry.
 
double getMeasure () const
 get measure of element
 
double & getMeasure ()
 get measure of element
 
MoFEM::Interface & getMField ()
 
moab::Interface & getMoab ()
 
virtual boost::weak_ptr< ForcesAndSourcesCore > getSubPipelinePtr () const
 
MoFEMErrorCode loopSide (const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t dim, const EntityHandle ent_for_side=0, boost::shared_ptr< Range > fe_range=nullptr, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy, AdjCache *adj_cache=nullptr)
 User calls this function to loop over elements on the side of face. This function calls finite element with its operator to do calculations.
 
MoFEMErrorCode loopThis (const string &fe_name, ForcesAndSourcesCore *this_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
 User calls this function to loop over the same element using a different set of integration points. This function calls finite element with its operator to do calculations.
 
MoFEMErrorCode loopParent (const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
 User calls this function to loop over parent elements. This function calls finite element with its operator to do calculations.
 
MoFEMErrorCode loopChildren (const string &fe_name, ForcesAndSourcesCore *child_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
 User calls this function to loop over parent elements. This function calls finite element with its operator to do calculations.
 
MoFEMErrorCode loopRange (const string &fe_name, ForcesAndSourcesCore *range_fe, boost::shared_ptr< Range > fe_range, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
 Iterate over range of elements.
 
- Public Member Functions inherited from MoFEM::DataOperator
 DataOperator (const bool symm=true)
 
virtual ~DataOperator ()=default
 
virtual MoFEMErrorCode doWork (int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
 Operator for bi-linear form, usually to calculate values on left hand side.
 
virtual MoFEMErrorCode opLhs (EntitiesFieldData &row_data, EntitiesFieldData &col_data)
 
virtual MoFEMErrorCode opRhs (EntitiesFieldData &data, const bool error_if_no_base=false)
 
bool getSymm () const
 Get if operator uses symmetry of DOFs or not.
 
void setSymm ()
 set if operator is executed taking in account symmetry
 
void unSetSymm ()
 unset if operator is executed for non symmetric problem
 

Private Attributes

boost::shared_ptr< DataAtIntegrationPts > dataAtPts
 data at integration pts
 

Additional Inherited Members

- Public Types inherited from MoFEM::ForcesAndSourcesCore::UserDataOperator
enum  OpType {
  OPROW = 1 << 0 , OPCOL = 1 << 1 , OPROWCOL = 1 << 2 , OPSPACE = 1 << 3 ,
  OPLAST = 1 << 3
}
 Controls loop over entities on element. More...
 
using AdjCache = std::map< EntityHandle, std::vector< boost::weak_ptr< NumeredEntFiniteElement > > >
 
- Public Types inherited from MoFEM::DataOperator
using DoWorkLhsHookFunType = boost::function< MoFEMErrorCode(DataOperator *op_ptr, int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)>
 
using DoWorkRhsHookFunType = boost::function< MoFEMErrorCode(DataOperator *op_ptr, int side, EntityType type, EntitiesFieldData::EntData &data)>
 
- Public Attributes inherited from MoFEM::ForcesAndSourcesCore::UserDataOperator
char opType
 
std::string rowFieldName
 
std::string colFieldName
 
FieldSpace sPace
 
- Public Attributes inherited from MoFEM::DataOperator
DoWorkLhsHookFunType doWorkLhsHook
 
DoWorkRhsHookFunType doWorkRhsHook
 
bool sYmm
 If true assume that matrix is symmetric structure.
 
std::array< bool, MBMAXTYPE > doEntities
 If true operator is executed for entity.
 
bool & doVertices
 \deprectaed If false skip vertices
 
bool & doEdges
 \deprectaed If false skip edges
 
bool & doQuads
 \deprectaed
 
bool & doTris
 \deprectaed
 
bool & doTets
 \deprectaed
 
bool & doPrisms
 \deprectaed
 
- Static Public Attributes inherited from MoFEM::ForcesAndSourcesCore::UserDataOperator
static const char *const OpTypeNames []
 
- Protected Member Functions inherited from MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
MoFEMErrorCode setPtrFE (ForcesAndSourcesCore *ptr)
 
- Protected Attributes inherited from MoFEM::ForcesAndSourcesCore::UserDataOperator
ForcesAndSourcesCore * ptrFE
 

Detailed Description

Examples
/home/lk58p/mofem_install/vanilla_dev_release/mofem-cephas/mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 287 of file EshelbianOperators.hpp.

Constructor & Destructor Documentation

◆ OpCalculateRotationAndSpatialGradient()

OpCalculateRotationAndSpatialGradient::OpCalculateRotationAndSpatialGradient ( boost::shared_ptr< DataAtIntegrationPts >  data_ptr)
inline

Definition at line 289 of file EshelbianOperators.hpp.

291 : VolUserDataOperator(NOSPACE, OPSPACE), dataAtPts(data_ptr) {}
NOSPACE
@ NOSPACE
Definition definitions.h:83
VolUserDataOperator
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
Definition elasticity.cpp:76
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPSPACE
@ OPSPACE
operator do Work is execute on space data
Definition ForcesAndSourcesCore.hpp:605
OpCalculateRotationAndSpatialGradient::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
Definition EshelbianOperators.hpp:296

Member Function Documentation

◆ doWork()

MoFEMErrorCode OpCalculateRotationAndSpatialGradient::doWork ( int  side,
EntityType  type,
EntData &  data 
)
virtual

Operator for linear form, usually to calculate values on right hand side.

Reimplemented from MoFEM::DataOperator.

Examples
/home/lk58p/mofem_install/vanilla_dev_release/mofem-cephas/mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianOperators.cpp.

Definition at line 175 of file EshelbianOperators.cpp.

177 {
178 MoFEMFunctionBegin;
179
180 auto ts_ctx = getTSCtx();
181 int nb_integration_pts = getGaussPts().size2();
182
183 // space size indices
184 FTENSOR_INDEX(SPACE_DIM, i);
185 FTENSOR_INDEX(SPACE_DIM, j);
186 FTENSOR_INDEX(SPACE_DIM, k);
187 FTENSOR_INDEX(SPACE_DIM, l);
188 FTENSOR_INDEX(SPACE_DIM, m);
189 FTENSOR_INDEX(SPACE_DIM, n);
190 FTENSOR_INDEX(SPACE_DIM, o);
191
192 // sym size indices
193 FTENSOR_INDEX(size_symm, L);
194
195 auto t_L = FTensor::SymmLTensor<double, 3>();
196
197 MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
198 *dataAtPts->getStretchTensorAtPts(), nb_integration_pts);
199 MatrixSizeHelper<GetFTensor4DdgFromMatType<3, 3, -1, DL>, DL>::size(
200 *dataAtPts->getDiffStretchTensorAtPts(), nb_integration_pts);
201 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
202 *dataAtPts->getStretchH1AtPts(), nb_integration_pts);
203 MatrixSizeHelper<GetFTensor4FromMatType<3, 3, 3, 3, -1, DL>, DL>::size(
204 *dataAtPts->getDiffStretchH1AtPts(), nb_integration_pts);
205 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
206 *dataAtPts->getAdjointPdstretchAtPts(), nb_integration_pts);
207 MatrixSizeHelper<GetFTensor1FromMatType<size_symm, -1, DL>, DL>::size(
208 *dataAtPts->getAdjointPdUAtPts(), nb_integration_pts);
209 MatrixSizeHelper<GetFTensor3FromMatType<3, 3, size_symm, -1, DL>, DL>::size(
210 *dataAtPts->getAdjointPdUdPAtPts(), nb_integration_pts);
211 MatrixSizeHelper<GetFTensor2FromMatType<3, size_symm, -1, DL>, DL>::size(
212 *dataAtPts->getAdjointPdUdOmegaAtPts(), nb_integration_pts);
213
214 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
215 *dataAtPts->getDeformationGradient(), nb_integration_pts);
216 MatrixSizeHelper<GetFTensor3FromMatType<3, 3, 3, -1, DL>, DL>::size(
217 dataAtPts->hdOmegaAtPts, nb_integration_pts);
218 MatrixSizeHelper<GetFTensor3FromMatType<3, 3, size_symm, -1, DL>, DL>::size(
219 dataAtPts->hdLogStretchAtPts, nb_integration_pts);
220
221 MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
222 dataAtPts->leviKirchhoffAtPts, nb_integration_pts);
223 MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
224 dataAtPts->leviKirchhoff0AtPts, nb_integration_pts);
225 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
226 dataAtPts->leviKirchhoffdOmegaAtPts, nb_integration_pts);
227 MatrixSizeHelper<GetFTensor2FromMatType<3, size_symm, -1, DL>, DL>::size(
228 dataAtPts->leviKirchhoffdLogStreatchAtPts, nb_integration_pts);
229 MatrixSizeHelper<GetFTensor3FromMatType<3, 3, 3, -1, DL>, DL>::size(
230 dataAtPts->leviKirchhoffPAtPts, nb_integration_pts);
231
232 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
233 dataAtPts->rotMatAtPts, nb_integration_pts);
234 MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
235 *dataAtPts->getEigenVals(), nb_integration_pts);
236 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
237 *dataAtPts->getEigenVecs(), nb_integration_pts);
238 dataAtPts->nbUniq.resize(nb_integration_pts, false);
239 MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>, DL>::size(
240 dataAtPts->eigenValsC, nb_integration_pts);
241 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
242 dataAtPts->eigenVecsC, nb_integration_pts);
243 dataAtPts->nbUniqC.resize(nb_integration_pts, false);
244
245 MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
246 dataAtPts->logStretch2H1AtPts, nb_integration_pts);
247 MatrixSizeHelper<GetFTensor2SymmetricFromMatType<3, -1, DL>, DL>::size(
248 dataAtPts->logStretchTotalTensorAtPts, nb_integration_pts);
249
250 MatrixSizeHelper<GetFTensor2FromMatType<3, 3, -1, DL>, DL>::size(
251 dataAtPts->internalStressAtPts, nb_integration_pts);
252 dataAtPts->internalStressAtPts.clear();
253
254 // Calculated values
255 auto t_h = dataAtPts->getFTensorSmallH(getGaussPts().size2());
256 auto t_h_domega = dataAtPts->getFTensorSmallHdOmega(getGaussPts().size2());
257 auto t_h_dlog_u =
258 dataAtPts->getFTensorSmallHdLogStretch(getGaussPts().size2());
259 auto t_levi_kirchhoff =
260 dataAtPts->getFTensorLeviKirchhoff(getGaussPts().size2());
261 auto t_levi_kirchhoff0 =
262 dataAtPts->getFTensorLeviKirchhoff0(getGaussPts().size2());
263 auto t_levi_kirchhoff_domega =
264 dataAtPts->getFTensorLeviKirchhoffdOmega(getGaussPts().size2());
265 auto t_levi_kirchhoff_dstreach =
266 dataAtPts->getFTensorLeviKirchhoffdLogStretch(getGaussPts().size2());
267 auto t_levi_kirchhoff_dP =
268 dataAtPts->getFTensorLeviKirchhoffP(getGaussPts().size2());
269 auto t_approx_P_adjoint_dstretch =
270 dataAtPts->getFTensorAdjointPdstretch(getGaussPts().size2());
271 auto t_approx_P_adjoint_log_du =
272 dataAtPts->getFTensorAdjointPdU(getGaussPts().size2());
273 auto t_approx_P_adjoint_log_du_dP =
274 dataAtPts->getFTensorAdjointPdUdP(getGaussPts().size2());
275 auto t_approx_P_adjoint_log_du_domega =
276 dataAtPts->getFTensorAdjointPdUdOmega(getGaussPts().size2());
277 auto t_R = dataAtPts->getFTensorRotMat(getGaussPts().size2());
278 auto t_u = dataAtPts->getFTensorStretch(getGaussPts().size2());
279 auto t_diff_u = dataAtPts->getFTensorDiffStretch(getGaussPts().size2());
280 auto t_eigen_vals = dataAtPts->getFTensorEigenVals(getGaussPts().size2());
281 auto t_eigen_vecs = dataAtPts->getFTensorEigenVecs(getGaussPts().size2());
282 auto &nbUniq = dataAtPts->nbUniq;
283 auto t_nb_uniq =
284 FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniq.data().data());
285 auto t_eigen_vals_C = dataAtPts->getFTensorEigenValsC(nb_integration_pts);
286 auto t_eigen_vecs_C = dataAtPts->getFTensorEigenVecsC(nb_integration_pts);
287 auto &nbUniqC = dataAtPts->nbUniqC;
288 auto t_nb_uniq_C =
289 FTensor::Tensor0<FTensor::PackPtr<int *, 1>>(nbUniqC.data().data());
290
291 auto t_u_h1 = dataAtPts->getFTensorStretchH1(getGaussPts().size2());
292 auto t_diff_u_h1 = dataAtPts->getFTensorDiffStretchH1(getGaussPts().size2());
293 auto t_log_stretch_total =
294 dataAtPts->getFTensorLogStretchTotal(getGaussPts().size2());
295 auto t_log_u2_h1 = dataAtPts->getFTensorLogStretch2H1(getGaussPts().size2());
296
297 // Field values
298 auto t_grad_h1 = dataAtPts->getFTensorSmallWGradH1(getGaussPts().size2());
299 auto t_omega = dataAtPts->getFTensorRotAxis(getGaussPts().size2());
300 auto t_approx_P = dataAtPts->getFTensorApproxP(getGaussPts().size2());
301 auto t_approx_P0 = dataAtPts->getFTensorApproxP0(getGaussPts().size2());
302 auto t_log_u = dataAtPts->getFTensorLogStretch(getGaussPts().size2());
303
304 // Rot axis 0
305 auto t_omega0 = dataAtPts->getFTensorRotAxis0(getGaussPts().size2());
306 auto t_log_u0 = dataAtPts->getFTensorLogStretch0(getGaussPts().size2());
307
308 auto next = [&]() {
309 // calculated values
310 ++t_h;
311 ++t_h_domega;
312 ++t_h_dlog_u;
313 ++t_levi_kirchhoff;
314 ++t_levi_kirchhoff0;
315 ++t_levi_kirchhoff_domega;
316 ++t_levi_kirchhoff_dstreach;
317 ++t_levi_kirchhoff_dP;
318 ++t_approx_P_adjoint_dstretch;
319 ++t_approx_P_adjoint_log_du;
320 ++t_approx_P_adjoint_log_du_dP;
321 ++t_approx_P_adjoint_log_du_domega;
322 ++t_R;
323 ++t_u;
324 ++t_diff_u;
325 ++t_eigen_vals;
326 ++t_eigen_vecs;
327 ++t_nb_uniq;
328 ++t_eigen_vals_C;
329 ++t_eigen_vecs_C;
330 ++t_nb_uniq_C;
331 ++t_u_h1;
332 ++t_diff_u_h1;
333 ++t_log_u2_h1;
334 ++t_log_stretch_total;
335 // field values
336 ++t_omega;
337 ++t_omega0;
338 ++t_grad_h1;
339 ++t_approx_P;
340 ++t_approx_P0;
341 ++t_log_u;
342 ++t_log_u0;
343 };
344
345 auto t_kd = FTensor::Kronecker_Delta<double>();
346 auto t_symm_kd = FTensor::Kronecker_Delta_symmetric<double>();
347 constexpr auto t_diff_sym = FTensor::DiffSymmetrize<double>();
348
349 auto calculate_stretch_from_log = [&](auto &t_log_u_src, auto &t_u_dst,
350 auto &t_eigen_vals_dst,
351 auto &t_eigen_vecs_dst,
352 int &nb_uniq_dst) {
353 MoFEMFunctionBegin;
354 FTensor::Tensor1<double, 3> eig;
355 FTensor::Tensor2<double, 3, 3> eigen_vec;
356 eigen_vec(i, j) = t_log_u_src(i, j);
357 if (computeEigenValuesSymmetric(eigen_vec, eig) != MB_SUCCESS) {
358 MOFEM_LOG("SELF", Sev::error) << "Failed to compute eigen values";
359 }
360 // CHKERR bound_eig(eig);
361 // rare case when two eigen values are equal
362 nb_uniq_dst = getUniqNb<3>(eig);
363 if (nb_uniq_dst < 3) {
364 CHKERR sortEigenVals<3>(eig, eigen_vec);
365 }
366 t_eigen_vals_dst(i) = eig(i);
367 t_eigen_vecs_dst(i, j) = eigen_vec(i, j);
368 t_u_dst(i, j) = EigenMatrix::getMat(t_eigen_vals_dst, t_eigen_vecs_dst,
369 EshelbianCore::f)(i, j);
370 MoFEMFunctionReturn(0);
371 };
372
373 auto calculate_log_stretch = [&]() {
374 MoFEMFunctionBegin;
375 int nb_uniq_val = 0;
376 CHKERR calculate_stretch_from_log(t_log_u, t_u, t_eigen_vals, t_eigen_vecs,
377 nb_uniq_val);
378 t_nb_uniq = nb_uniq_val;
379 auto get_t_diff_u = [&]() {
380 return EigenMatrix::getDiffMat(t_eigen_vals, t_eigen_vecs,
381 EshelbianCore::f, EshelbianCore::d_f,
382 t_nb_uniq);
383 };
384 t_diff_u(i, j, k, l) = get_t_diff_u()(i, j, k, l);
385 FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
386 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
387 MoFEMFunctionReturn(0);
388 };
389
390 auto calculate_total_stretch = [&](auto &t_h1) {
391 MoFEMFunctionBegin;
392 if (EshelbianCore::gradApproximator == NO_H1_CONFIGURATION) {
393
394 t_log_u2_h1(i, j) = 0;
395 t_log_stretch_total(i, j) = t_log_u(i, j);
396
397 } else {
398
399 FTensor::Tensor1<double, 3> t_eig_C;
400 FTensor::Tensor1<double, 3> t_coordinate_stretch;
401 FTensor::Tensor2<double, 3, 3> t_eigen_vec;
402
403 FTensor::Tensor2<double, 3, 3> t_C_h1;
404 t_C_h1(i, j) = t_h1(k, i) * t_h1(k, j);
405 t_eigen_vec(i, j) = t_C_h1(i, j);
406 if (computeEigenValuesSymmetric(t_eigen_vec, t_eig_C) != MB_SUCCESS) {
407 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
408 "Failed to compute eigenvalues of F_H1^T F_H1");
409 }
410 // rare case when two eigen values are equal
411 t_nb_uniq_C = getUniqNb<3>(t_eig_C);
412 if (t_nb_uniq_C < 3) {
413 CHKERR sortEigenVals<3>(t_eig_C, t_eigen_vec);
414 }
415 for (int aa = 0; aa != 3; ++aa) {
416 if (!std::isfinite(t_eig_C(aa)) || t_eig_C(aa) <= 0.) {
417 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
418 "F_H1^T F_H1 must be positive definite; eigenvalue %d is "
419 "%g",
420 aa, t_eig_C(aa));
421 }
422 const double principal_stretch = std::sqrt(t_eig_C(aa));
423 const double coordinate_stretch =
424 EshelbianCore::inv_f(principal_stretch);
425 if (!std::isfinite(coordinate_stretch)) {
426 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP,
427 "Non-finite H1 coordinate stretch for principal stretch %g",
428 principal_stretch);
429 }
430 t_coordinate_stretch(aa) = coordinate_stretch;
431 }
432 t_eigen_vals_C(i) = t_eig_C(i);
433 t_eigen_vecs_C(i, j) = t_eigen_vec(i, j);
434
435 t_log_u2_h1(i, j) =
436 EigenMatrix::getMat(t_coordinate_stretch, t_eigen_vec,
437 [](const double v) { return v; })(i, j);
438 // The hand-coded Hencky formulation uses additive stretch coordinates.
439 // For logarithmic coordinates this is log(U_H1) + log(U_increment).
440 t_log_stretch_total(i, j) = t_log_u2_h1(i, j) + t_log_u(i, j);
441 }
442 MoFEMFunctionReturn(0);
443 };
444
445 auto no_h1_loop = [&]() {
446 MoFEMFunctionBegin;
447
448 switch (EshelbianCore::rotSelector) {
449 case LARGE_ROT:
450 break;
451 case SMALL_ROT:
452 break;
453 default:
454 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
455 "no_h1_loop is only implemented for LARGE_ROT");
456 };
457
458 for (int gg = 0; gg != nb_integration_pts; ++gg) {
459
460 auto t_kd = FTensor::Kronecker_Delta<double>();
461
462 FTensor::Tensor2<double, 3, 3> t_h1;
463 t_h1(i, j) = t_kd(i, j);
464
465 // calculate streach
466 CHKERR calculate_log_stretch();
467 FTensor::Tensor2<double, 3, 3> t_u0;
468 if (EshelbianCore::isNoStretch()) {
469 t_u0(i, j) = t_u(i, j);
470 } else {
471 FTensor::Tensor1<double, 3> t_eigen_vals_0;
472 FTensor::Tensor2<double, 3, 3> t_eigen_vecs_0;
473 int nb_uniq_0 = 0;
474 CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
475 t_eigen_vecs_0, nb_uniq_0);
476 }
477 // calculate total stretch
478 CHKERR calculate_total_stretch(t_h1);
479
480 t_u_h1(i, j) = t_u(i, j);
481 t_diff_u_h1(i, j, k, l) = t_diff_u(i, j, k, l);
482 FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
483 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
484
485 FTensor::Tensor3<double, 3, 3, 3> t_diff_R;
486 FTensor::Tensor4<double, 3, 3, 3, 3> t_diff_diff_R;
487
488 auto large_rot = [&]() {
489 t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
490 t_diff_R(i, j, k) =
491 LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
492 t_diff_diff_R(i, j, k, l) =
493 LieGroups::SO3::diffDiffExp(t_omega, t_omega.l2())(i, j, k, l);
494
495 FTensor::Tensor3<double, 3, 3, 3> t_diff_R0;
496 t_diff_R0(i, j, k) =
497 LieGroups::SO3::diffExp(t_omega0, t_omega0.l2())(i, j, k);
498
499 t_h(i, k) = t_R(i, l) * t_u(l, k);
500
501 FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
502 t_rotated_P(l, k) = t_R(i, l) * t_approx_P(i, k);
503 t_approx_P_adjoint_dstretch(l, k) =
504 t_diff_sym(l, k, i, j) * t_rotated_P(i, j);
505 t_approx_P_adjoint_log_du(L) =
506 t_approx_P_adjoint_dstretch(l, k) * t_Ldiff_u(l, k, L);
507
508 t_levi_kirchhoff(m) =
509 t_diff_R(i, l, m) * (t_u(l, k) * t_approx_P(i, k));
510 t_levi_kirchhoff0(m) =
511 t_diff_R0(i, l, m) * (t_u0(l, k) * t_approx_P0(i, k));
512
513 if (ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
514 t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_u(l, k);
515 t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_u(l, k, L);
516
517 t_approx_P_adjoint_log_du_dP(i, k, L) =
518 t_R(i, l) * t_Ldiff_u(l, k, L);
519
520 FTensor::Tensor3<double, 3, 3, 3> t_A;
521 t_A(k, l, m) = t_diff_R(i, l, m) * t_approx_P(i, k);
522 t_approx_P_adjoint_log_du_domega(m, L) =
523 t_A(k, l, m) * t_Ldiff_u(k, l, L);
524
525 t_levi_kirchhoff_dstreach(m, L) =
526 t_diff_R(i, l, m) * (t_Ldiff_u(l, k, L) * t_approx_P(i, k));
527 t_levi_kirchhoff_dP(m, i, k) = t_diff_R(i, l, m) * t_u(l, k);
528 t_levi_kirchhoff_domega(m, n) =
529 t_diff_diff_R(i, l, m, n) * (t_u(l, k) * t_approx_P(i, k));
530
531 if (EshelbianCore::stretchHandling ==
532 EshelbianCore::NO_STREACH_NONLINEAR) {
533 auto t_d_u_d_b = GetFTensor4DdgFromMatImpl<
534 SPACE_DIM, SPACE_DIM, -1, DL, MatrixDouble>::get(
535 dataAtPts->matInvD, gg, 0);
536 auto [t_d_b_d_omega, t_d_u_d_omega, t_d_h_d_omega] =
537 getDiffSpatialGradientDR(t_d_u_d_b, t_R, t_diff_R,
538 t_approx_P);
539
540 t_h_domega(i, k, m) += t_d_h_d_omega(i, k, m);
541 FTensor::Tensor3<double, SPACE_DIM, SPACE_DIM, SPACE_DIM>
542 t_d_u_contract_p;
543 t_d_u_contract_p(i, l, n) =
544 t_d_u_d_omega(l, k, n) * t_approx_P(i, k);
545 t_levi_kirchhoff_domega(m, n) +=
546 t_diff_R(i, l, m) * t_d_u_contract_p(i, l, n);
547
548 if constexpr (EshelbianCore::symmetrySelector > SYMMETRIC) {
549 FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM,
550 SPACE_DIM>
551 t_d_b_d_p;
552 t_d_b_d_p(i, j, k, l) =
553 t_diff_sym(i, j, m, l) * t_R(k, m);
554 FTensor::Tensor4<double, SPACE_DIM, SPACE_DIM, SPACE_DIM,
555 SPACE_DIM>
556 t_d_u_d_p;
557 t_d_u_d_p(i, j, k, l) =
558 t_d_u_d_b(i, j, m, n) * t_d_b_d_p(m, n, k, l);
559 t_levi_kirchhoff_dP(m, k, l) +=
560 t_d_u_d_p(i, j, k, l) * t_d_b_d_omega(i, j, m);
561 }
562 }
563 }
564 };
565
566 auto moderate_rot = [&](auto &t_omega0) {
567 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
568 "moderate_rot is not implemented yet");
569 };
570
571 auto small_rot = [&]() {
572 t_u_h1(i, j) = t_u(i, j);
573 t_diff_u_h1(i, j, k, l) = t_diff_u(i, j, k, l);
574 FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
575 t_Ldiff_u(i, j, L) = t_diff_u(i, j, m, n) * t_L(m, n, L);
576
577 t_R(i, j) = t_kd(i, j) + levi_civita(i, j, k) * t_omega(k);
578 t_h(i, j) = levi_civita(i, j, k) * t_omega(k) + t_u(i, j);
579
580 t_h_domega(i, j, k) = levi_civita(i, j, k);
581 t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
582
583 // Adjoint stress
584 FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
585 t_rotated_P(i, j) = t_R(k, i) * t_approx_P(k, j);
586 t_approx_P_adjoint_dstretch(i, j) =
587 t_diff_sym(i, j, k, l) * t_rotated_P(k, l);
588 t_approx_P_adjoint_log_du(L) =
589 t_approx_P_adjoint_dstretch(i, j) * t_Ldiff_u(i, j, L);
590 t_approx_P_adjoint_log_du_dP(i, j, L) = t_Ldiff_u(i, j, L);
591 t_approx_P_adjoint_log_du_domega(m, L) = 0;
592
593 // Kirchhoff stress
594 t_levi_kirchhoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
595 t_levi_kirchhoff0(k) = levi_civita(i, j, k) * t_approx_P0(i, j);
596 t_levi_kirchhoff_dstreach(m, L) = 0;
597 t_levi_kirchhoff_dP(k, i, j) = levi_civita(i, j, k);
598 t_levi_kirchhoff_domega(m, n) = 0;
599 };
600
601 // rotation
602 switch (EshelbianCore::rotSelector) {
603 case LARGE_ROT:
604 large_rot();
605 break;
606 case MODERATE_ROT:
607 moderate_rot(t_omega0);
608 break;
609 case SMALL_ROT:
610 small_rot();
611 break;
612 default:
613 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
614 "rotationSelector not handled");
615 }
616
617 next();
618 }
619
620 MoFEMFunctionReturn(0);
621 };
622
623 auto large_loop = [&]() {
624 MoFEMFunctionBegin;
625
626 switch (EshelbianCore::rotSelector) {
627 case LARGE_ROT:
628 break;
629 case SMALL_ROT:
630 break;
631 default:
632 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
633 "rotSelector should be large or small");
634 };
635
636 for (int gg = 0; gg != nb_integration_pts; ++gg) {
637
638 auto t_kd = FTensor::Kronecker_Delta<double>();
639
640 FTensor::Tensor2<double, 3, 3> t_h1;
641 switch (EshelbianCore::gradApproximator) {
642 case LARGE_ROT:
643 t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
644 break;
645 default:
646 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
647 "Selected grad approximator not handled");
648 };
649
650 // calculate streach
651 CHKERR calculate_log_stretch();
652 FTensor::Tensor2<double, 3, 3> t_u0;
653 if (EshelbianCore::isNoStretch()) {
654 t_u0(i, j) = t_u(i, j);
655 } else {
656 FTensor::Tensor1<double, 3> t_eigen_vals_0;
657 FTensor::Tensor2<double, 3, 3> t_eigen_vecs_0;
658 int nb_uniq_0 = 0;
659 CHKERR calculate_stretch_from_log(t_log_u0, t_u0, t_eigen_vals_0,
660 t_eigen_vecs_0, nb_uniq_0);
661 }
662 // calculate total stretch
663 CHKERR calculate_total_stretch(t_h1);
664
665 t_u_h1(l, k) = t_u(l, o) * t_h1(o, k);
666 FTensor::Tensor2<double, 3, 3> t_u_h10;
667 t_u_h10(l, k) = t_u0(l, o) * t_h1(o, k);
668 t_diff_u_h1(i, j, k, l) = t_diff_u(i, o, k, l) * t_h1(o, j);
669 FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u_h1;
670 t_Ldiff_u_h1(l, k, L) = t_diff_u_h1(l, k, i, j) * t_L(i, j, L);
671
672 FTensor::Tensor3<double, 3, 3, 3> t_diff_R;
673 FTensor::Tensor3<double, 3, 3, 3> t_diff_R0;
674 FTensor::Tensor4<double, 3, 3, 3, 3> t_diff_diff_R;
675
676 // rotation
677 switch (EshelbianCore::rotSelector) {
678 case SMALL_ROT:
679 t_R(i, k) = t_kd(i, k) + levi_civita(i, k, l) * t_omega(l);
680 t_diff_R(i, j, k) = levi_civita(i, j, k);
681 t_diff_R0(i, j, k) = levi_civita(i, j, k);
682 t_diff_diff_R(i, j, l, m) = 0;
683 break;
684 case LARGE_ROT:
685 t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
686 t_diff_R(i, j, k) =
687 LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
688 t_diff_R0(i, j, k) =
689 LieGroups::SO3::diffExp(t_omega0, t_omega0.l2())(i, j, k);
690 t_diff_diff_R(i, j, k, l) =
691 LieGroups::SO3::diffDiffExp(t_omega, t_omega.l2())(i, j, k, l);
692 break;
693
694 default:
695 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
696 "rotationSelector not handled");
697 }
698
699 // calculate gradient
700 t_h(i, k) = t_R(i, l) * t_u_h1(l, k);
701
702 // Adjoint stress
703 FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
704 t_rotated_P(l, o) =
705 (t_R(i, l) * t_approx_P(i, k)) * t_h1(o, k);
706 t_approx_P_adjoint_dstretch(l, o) =
707 t_diff_sym(l, o, i, j) * t_rotated_P(i, j);
708 t_approx_P_adjoint_log_du(L) =
709 t_R(i, l) * t_approx_P(i, k) * t_Ldiff_u_h1(l, k, L);
710
711 // Kirchhoff stress
712 t_levi_kirchhoff(m) = t_diff_R(i, l, m) * t_u_h1(l, k) * t_approx_P(i, k);
713 t_levi_kirchhoff0(m) =
714 t_diff_R0(i, l, m) * t_u_h10(l, k) * t_approx_P0(i, k);
715
716 if (ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
717
718 t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_u_h1(l, k);
719 t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_u_h1(l, k, L);
720
721 t_approx_P_adjoint_log_du_dP(i, k, L) =
722 t_R(i, l) * t_Ldiff_u_h1(l, k, L);
723
724 FTensor::Tensor4<double, 3, size_symm, 3, 3> t_A;
725 t_A(m, L, i, k) = t_diff_R(i, l, m) * t_Ldiff_u_h1(l, k, L);
726 t_approx_P_adjoint_log_du_domega(m, L) =
727 t_A(m, L, i, k) * t_approx_P(i, k);
728
729 t_levi_kirchhoff_dstreach(m, L) =
730 t_diff_R(i, l, m) * (t_Ldiff_u_h1(l, k, L) * t_approx_P(i, k));
731
732 t_levi_kirchhoff_dP(m, i, k) = t_diff_R(i, l, m) * t_u_h1(l, k);
733 t_levi_kirchhoff_domega(m, n) =
734 t_diff_diff_R(i, l, m, n) * (t_u_h1(l, k) * t_approx_P(i, k));
735 }
736
737 next();
738 }
739
740 MoFEMFunctionReturn(0);
741 };
742
743 auto moderate_loop = [&]() {
744 MoFEMFunctionBegin;
745
746 switch (EshelbianCore::rotSelector) {
747 case LARGE_ROT:
748 break;
749 case SMALL_ROT:
750 break;
751 default:
752 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
753 "rotSelector should be large or small");
754 };
755
756 for (int gg = 0; gg != nb_integration_pts; ++gg) {
757
758 auto t_kd = FTensor::Kronecker_Delta<double>();
759
760 FTensor::Tensor2<double, 3, 3> t_h1;
761 switch (EshelbianCore::gradApproximator) {
762 case MODERATE_ROT:
763 t_h1(i, j) = t_grad_h1(i, j) + t_kd(i, j);
764 break;
765 default:
766 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
767 "Selected grad approximator not handled");
768 };
769
770 // calculate streach
771 CHKERR calculate_log_stretch();
772 // calculate total stretch
773 CHKERR calculate_total_stretch(t_h1);
774
775 auto t_diff = FTensor::DiffTensor<double>();
776
777 t_u_h1(l, k) = (t_kd(l, o) + t_log_u(l, o)) * t_h1(o, k);
778 FTensor::Tensor2<double, 3, 3> t_u_h10;
779 t_u_h10(l, k) = (t_kd(l, o) + t_log_u0(l, o)) * t_h1(o, k);
780 t_diff_u_h1(i, j, k, l) = t_diff(i, o, k, l) * t_h1(o, j);
781 FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u_h1;
782 t_Ldiff_u_h1(l, k, L) = t_diff_u_h1(l, k, i, j) * t_L(i, j, L);
783
784 FTensor::Tensor3<double, 3, 3, 3> t_diff_R;
785 FTensor::Tensor3<double, 3, 3, 3> t_diff_R0;
786 FTensor::Tensor4<double, 3, 3, 3, 3> t_diff_diff_R;
787
788 // rotation
789 switch (EshelbianCore::rotSelector) {
790 case SMALL_ROT:
791 t_R(i, k) = t_kd(i, k) + levi_civita(i, k, l) * t_omega(l);
792 t_diff_R(i, j, k) = levi_civita(i, j, k);
793 t_diff_R0(i, j, k) = levi_civita(i, j, k);
794 t_diff_diff_R(i, j, l, m) = 0;
795 break;
796 case LARGE_ROT:
797 t_R(i, j) = LieGroups::SO3::exp(t_omega, t_omega.l2())(i, j);
798 t_diff_R(i, j, k) =
799 LieGroups::SO3::diffExp(t_omega, t_omega.l2())(i, j, k);
800 t_diff_R0(i, j, k) =
801 LieGroups::SO3::diffExp(t_omega0, t_omega0.l2())(i, j, k);
802 t_diff_diff_R(i, j, k, l) =
803 LieGroups::SO3::diffDiffExp(t_omega, t_omega.l2())(i, j, k, l);
804 break;
805
806 default:
807 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
808 "rotationSelector not handled");
809 }
810
811 // calculate gradient
812 t_h(i, k) = t_R(i, l) * t_u_h1(l, k);
813
814 // Adjoint stress
815 FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_rotated_P;
816 t_rotated_P(l, o) =
817 (t_R(i, l) * t_approx_P(i, k)) * t_h1(o, k);
818 t_approx_P_adjoint_dstretch(l, o) =
819 t_diff_sym(l, o, i, j) * t_rotated_P(i, j);
820 t_approx_P_adjoint_log_du(L) =
821 t_R(i, l) * t_approx_P(i, k) * t_Ldiff_u_h1(l, k, L);
822
823 // Kirchhoff stress
824 t_levi_kirchhoff(m) = t_diff_R(i, l, m) * t_u_h1(l, k) * t_approx_P(i, k);
825 t_levi_kirchhoff0(m) =
826 t_diff_R0(i, l, m) * t_u_h10(l, k) * t_approx_P0(i, k);
827
828 if (ts_ctx == TSMethod::CTX_TSSETIJACOBIAN) {
829
830 t_h_domega(i, k, m) = t_diff_R(i, l, m) * t_u_h1(l, k);
831 t_h_dlog_u(i, k, L) = t_R(i, l) * t_Ldiff_u_h1(l, k, L);
832
833 t_approx_P_adjoint_log_du_dP(i, k, L) =
834 t_R(i, l) * t_Ldiff_u_h1(l, k, L);
835
836 FTensor::Tensor4<double, 3, size_symm, 3, 3> t_A;
837 t_A(m, L, i, k) = t_diff_R(i, l, m) * t_Ldiff_u_h1(l, k, L);
838 t_approx_P_adjoint_log_du_domega(m, L) =
839 t_A(m, L, i, k) * t_approx_P(i, k);
840
841 t_levi_kirchhoff_dstreach(m, L) =
842 t_diff_R(i, l, m) * (t_Ldiff_u_h1(l, k, L) * t_approx_P(i, k));
843
844 t_levi_kirchhoff_dP(m, i, k) = t_diff_R(i, l, m) * t_u_h1(l, k);
845 t_levi_kirchhoff_domega(m, n) =
846 t_diff_diff_R(i, l, m, n) * (t_u_h1(l, k) * t_approx_P(i, k));
847 }
848
849 next();
850 }
851
852 MoFEMFunctionReturn(0);
853 };
854
855 auto small_loop = [&]() {
856 MoFEMFunctionBegin;
857 switch (EshelbianCore::rotSelector) {
858 case SMALL_ROT:
859 break;
860 default:
861 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
862 "rotSelector should be small");
863 };
864
865 for (int gg = 0; gg != nb_integration_pts; ++gg) {
866
867 FTensor::Tensor2<double, 3, 3> t_h1;
868 switch (EshelbianCore::gradApproximator) {
869 case SMALL_ROT:
870 t_h1(i, j) = t_kd(i, j);
871 break;
872 default:
873 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
874 "gradApproximator not handled");
875 };
876
877 FTensor::Tensor3<double, 3, 3, size_symm> t_Ldiff_u;
878 if (EshelbianCore::stretchSelector > LINEAR) {
879 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
880 "stretchSelector should be linear for small loop");
881 } else {
882 t_u(i, j) = t_symm_kd(i, j) + t_log_u(i, j);
883 t_u_h1(i, j) = t_u(i, j);
884 t_diff_u_h1(i, j, k, l) =
885 (t_kd(i, k) * t_kd(j, l) + t_kd(i, l) * t_kd(j, k));
886 t_diff_u_h1(i, j, k, l) /= 2.;
887 t_Ldiff_u(i, j, L) = t_L(i, j, L);
888 }
889 t_log_u2_h1(i, j) = 0;
890 t_log_stretch_total(i, j) = t_log_u(i, j);
891
892 t_R(i, j) = t_kd(i, j) + levi_civita(i, j, k) * t_omega(k);
893 t_h(i, j) = levi_civita(i, j, k) * t_omega(k) + t_u(i, j);
894
895 t_h_domega(i, j, k) = levi_civita(i, j, k);
896 t_h_dlog_u(i, j, L) = t_Ldiff_u(i, j, L);
897
898 // Adjoint stress
899 t_approx_P_adjoint_dstretch(i, j) =
900 t_diff_sym(i, j, k, l) * t_approx_P(k, l);
901 t_approx_P_adjoint_log_du(L) =
902 t_approx_P_adjoint_dstretch(i, j) * t_Ldiff_u(i, j, L);
903 t_approx_P_adjoint_log_du_dP(i, j, L) = t_Ldiff_u(i, j, L);
904 t_approx_P_adjoint_log_du_domega(m, L) = 0;
905
906 // Kirchhoff stress
907 t_levi_kirchhoff(k) = levi_civita(i, j, k) * t_approx_P(i, j);
908 t_levi_kirchhoff0(k) = levi_civita(i, j, k) * t_approx_P0(i, j);
909 t_levi_kirchhoff_dstreach(m, L) = 0;
910 t_levi_kirchhoff_dP(k, i, j) = levi_civita(i, j, k);
911 t_levi_kirchhoff_domega(m, n) = 0;
912
913 next();
914 }
915
916 MoFEMFunctionReturn(0);
917 };
918
919 switch (EshelbianCore::gradApproximator) {
920 case NO_H1_CONFIGURATION:
921 CHKERR no_h1_loop();
922 break;
923 case LARGE_ROT:
924 CHKERR large_loop();
925 MoFEMFunctionReturnHot(0);
926 break;
927 case MODERATE_ROT:
928 CHKERR moderate_loop();
929 MoFEMFunctionReturnHot(0);
930 break;
931 case SMALL_ROT:
932 CHKERR small_loop();
933 MoFEMFunctionReturnHot(0);
934 break;
935 default:
936 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
937 "gradApproximator not handled");
938 break;
939 };
940
941 MoFEMFunctionReturn(0);
942}
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2673
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
FTensor::DiffSymmetrize
Fourth-order symmetrization tensor.
Definition DiffSymmetrize.hpp:15
FTensor::DiffTensor
Fourth-order differential tensor symmetric in both index pairs.
Definition DiffTensor.hpp:15
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition Kronecker_Delta.hpp:49
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15
FTensor::SymmLTensor
Mapping from symmetric tensor indices to packed storage index.
Definition SymmLTensor.hpp:16
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
t_kd
constexpr auto t_kd
Definition free_surface.cpp:161
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:41
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:39
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition level_set.cpp:1923
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:23
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
EigenMatrix::getMat
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
Definition MatrixFunction.hpp:133
EigenMatrix::getDiffMat
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition MatrixFunction.hpp:201
EshelbianPlasticity::getDiffSpatialGradientDR
auto getDiffSpatialGradientDR(TInvD &t_d_u_d_b, TRotation &t_R, TDiffRotation &t_diff_R, TStress &t_P)
Definition EshelbianOperators.cpp:155
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition Levi_Civita.hpp:617
MatOps::DL
DataLayoutTraits< DataLayout::GaussByCoeffs > DL
Definition MatHuHu.hpp:33
MoFEM::GetFTensor2SymmetricFromMatType
decltype(GetFTensor2SymmetricFromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2SymmetricFromMatType
Definition Templates.hpp:836
MoFEM::GetFTensor4FromMatType
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
Definition Templates.hpp:1360
MoFEM::GetFTensor4DdgFromMatType
decltype(GetFTensor4DdgFromMatImpl< Tensor_Dim01, Tensor_Dim23, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4DdgFromMatType
Definition Templates.hpp:1540
MoFEM::GetFTensor1FromMatType
decltype(GetFTensor1FromMatImpl< Tensor_Dim, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor1FromMatType
Definition Templates.hpp:385
MoFEM::computeEigenValuesSymmetric
MoFEMErrorCode computeEigenValuesSymmetric(const MatrixDouble &mat, VectorDouble &eig, MatrixDouble &eigen_vec)
compute eigenvalues of a symmetric matrix using lapack dsyev
Definition Templates.hpp:2118
MoFEM::GetFTensor3FromMatType
decltype(GetFTensor3FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor3FromMatType
Definition Templates.hpp:1013
MoFEM::GetFTensor2FromMatType
decltype(GetFTensor2FromMatImpl< Tensor_Dim0, Tensor_Dim1, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor2FromMatType
Definition Templates.hpp:551
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
EshelbianCore::stretchSelector
static enum StretchSelector stretchSelector
Definition EshelbianCore.hpp:54
EshelbianCore::stretchHandling
static enum StretchHandling stretchHandling
Definition EshelbianCore.hpp:55
EshelbianCore::rotSelector
static enum RotSelector rotSelector
Definition EshelbianCore.hpp:52
EshelbianCore::gradApproximator
static enum RotSelector gradApproximator
Definition EshelbianCore.hpp:53
EshelbianCore::symmetrySelector
static constexpr enum SymmetrySelector symmetrySelector
Definition EshelbianCore.hpp:48
EshelbianCore::f
static boost::function< double(const double)> f
Definition EshelbianCore.hpp:131
EshelbianCore::d_f
static boost::function< double(const double)> d_f
Definition EshelbianCore.hpp:132
EshelbianCore::isNoStretch
static bool isNoStretch()
Definition EshelbianCore.hpp:57
EshelbianCore::inv_f
static boost::function< double(const double)> inv_f
Definition EshelbianCore.hpp:134
LieGroups::SO3::diffDiffExp
static auto diffDiffExp(A &&t_w_vee, B &&theta)
Definition Lie.hpp:110
LieGroups::SO3::diffExp
static auto diffExp(A &&t_w_vee, B &&theta)
Definition Lie.hpp:105
LieGroups::SO3::exp
static auto exp(A &&t_w_vee, B &&theta)
Definition Lie.hpp:69
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition ForcesAndSourcesCore.hpp:1350
MoFEM::TSMethod::CTX_TSSETIJACOBIAN
@ CTX_TSSETIJACOBIAN
Setting up implicit Jacobian.
Definition LoopMethods.hpp:275
size_symm
constexpr auto size_symm
Definition plastic.cpp:42

Member Data Documentation

◆ dataAtPts

boost::shared_ptr<DataAtIntegrationPts> OpCalculateRotationAndSpatialGradient::dataAtPts
private

data at integration pts

Definition at line 296 of file EshelbianOperators.hpp.

The documentation for this struct was generated from the following files:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow