13using namespace boost::numeric;
25 boost::shared_ptr<MatrixDouble> mat_coords_ptr,
26 boost::shared_ptr<VectorDouble> density_at_pts,
27 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr)
33 HookeElement::EntData &row_data) {
35 if (row_type != MBVERTEX)
38 const int nb_integration_pts = getGaussPts().size2();
51 coords = trans(getCoordsAtGaussPts());
54 auto t_coords = getFTensor1FromMat<3>(coords);
56 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
57 t_rho = 1 + t_coords(
i) * t_coords(
i);
59 t_grad_rho(
i) = 2 * t_coords(
i);
70int main(
int argc,
char *argv[]) {
76 moab::Core mb_instance;
77 moab::Interface &moab = mb_instance;
84 PetscBool ale = PETSC_FALSE;
86 PetscBool test_jacobian = PETSC_FALSE;
100 if (ale == PETSC_TRUE) {
118 if (ale == PETSC_TRUE) {
124 boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
126 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
131 fe_lhs_ptr->getRuleHook =
VolRule();
132 fe_rhs_ptr->getRuleHook =
VolRule();
140 boost::shared_ptr<map<int, BlockData>> block_sets_ptr =
141 boost::make_shared<map<int, BlockData>>();
142 (*block_sets_ptr)[0].
iD = 0;
143 (*block_sets_ptr)[0].E = 1;
144 (*block_sets_ptr)[0].PoissonRatio = 0.25;
149 const double rho_n = 2.0;
150 const double rho_0 = 0.5;
152 auto my_operators = [&](boost::shared_ptr<ForcesAndSourcesCore> &fe_lhs_ptr,
153 boost::shared_ptr<ForcesAndSourcesCore> &fe_rhs_ptr,
154 boost::shared_ptr<map<int, BlockData>>
156 const std::string x_field,
157 const std::string X_field,
const bool ale,
158 const bool field_disp) {
161 boost::shared_ptr<HookeElement::DataAtIntegrationPts> data_at_pts(
162 new HookeElement::DataAtIntegrationPts());
163 boost::shared_ptr<MatrixDouble> mat_coords_ptr =
164 boost::make_shared<MatrixDouble>();
165 boost::shared_ptr<VectorDouble> rho_at_gauss_pts_ptr =
166 boost::make_shared<VectorDouble>();
167 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr =
168 boost::make_shared<MatrixDouble>();
171 if (ale == PETSC_FALSE) {
172 fe_lhs_ptr->getOpPtrVector().push_back(
175 x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
176 rho_grad_at_gauss_pts_ptr));
177 fe_lhs_ptr->getOpPtrVector().push_back(
178 new HookeElement::OpCalculateStiffnessScaledByDensityField(
179 x_field, x_field, block_sets_ptr, data_at_pts,
180 rho_at_gauss_pts_ptr, rho_n, rho_0));
181 fe_lhs_ptr->getOpPtrVector().push_back(
182 new HookeElement::OpLhs_dx_dx<1>(x_field, x_field, data_at_pts));
184 fe_lhs_ptr->getOpPtrVector().push_back(
187 fe_lhs_ptr->getOpPtrVector().push_back(
190 X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
191 rho_grad_at_gauss_pts_ptr));
192 fe_lhs_ptr->getOpPtrVector().push_back(
193 new HookeElement::OpCalculateStiffnessScaledByDensityField(
194 x_field, x_field, block_sets_ptr, data_at_pts,
195 rho_at_gauss_pts_ptr, rho_n, rho_0));
196 fe_lhs_ptr->getOpPtrVector().push_back(
199 fe_lhs_ptr->getOpPtrVector().push_back(
200 new HookeElement::OpCalculateStrainAle(x_field, x_field,
202 fe_lhs_ptr->getOpPtrVector().push_back(
203 new HookeElement::OpCalculateStress<1>(x_field, x_field,
205 fe_lhs_ptr->getOpPtrVector().push_back(
206 new HookeElement::OpAleLhs_dx_dx<1>(x_field, x_field,
208 fe_lhs_ptr->getOpPtrVector().push_back(
209 new HookeElement::OpAleLhs_dx_dX<1>(x_field, X_field,
211 fe_lhs_ptr->getOpPtrVector().push_back(
212 new HookeElement::OpCalculateEnergy(X_field, X_field,
214 fe_lhs_ptr->getOpPtrVector().push_back(
215 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
217 fe_lhs_ptr->getOpPtrVector().push_back(
218 new HookeElement::OpAleLhs_dX_dX<1>(X_field, X_field,
220 fe_lhs_ptr->getOpPtrVector().push_back(
221 new HookeElement::OpAleLhsPre_dX_dx<1>(X_field, x_field,
223 fe_lhs_ptr->getOpPtrVector().push_back(
224 new HookeElement::OpAleLhs_dX_dx(X_field, x_field, data_at_pts));
225 fe_lhs_ptr->getOpPtrVector().push_back(
226 new HookeElement::OpAleLhsWithDensity_dx_dX(
227 x_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
228 rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
229 fe_lhs_ptr->getOpPtrVector().push_back(
230 new HookeElement::OpAleLhsWithDensity_dX_dX(
231 X_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
232 rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
238 if (ale == PETSC_FALSE) {
239 fe_rhs_ptr->getOpPtrVector().push_back(
242 fe_rhs_ptr->getOpPtrVector().push_back(
245 x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
246 rho_grad_at_gauss_pts_ptr));
247 fe_rhs_ptr->getOpPtrVector().push_back(
248 new HookeElement::OpCalculateStiffnessScaledByDensityField(
249 x_field, x_field, block_sets_ptr, data_at_pts,
250 rho_at_gauss_pts_ptr, rho_n, rho_0));
252 fe_rhs_ptr->getOpPtrVector().push_back(
253 new HookeElement::OpCalculateStrain<1>(x_field, x_field,
256 fe_rhs_ptr->getOpPtrVector().push_back(
257 new HookeElement::OpCalculateStrain<0>(x_field, x_field,
260 fe_rhs_ptr->getOpPtrVector().push_back(
261 new HookeElement::OpCalculateStress<1>(x_field, x_field,
263 fe_rhs_ptr->getOpPtrVector().push_back(
264 new HookeElement::OpRhs_dx(x_field, x_field, data_at_pts));
266 fe_rhs_ptr->getOpPtrVector().push_back(
269 fe_rhs_ptr->getOpPtrVector().push_back(
272 X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
273 rho_grad_at_gauss_pts_ptr));
274 fe_rhs_ptr->getOpPtrVector().push_back(
275 new HookeElement::OpCalculateStiffnessScaledByDensityField(
276 x_field, x_field, block_sets_ptr, data_at_pts,
277 rho_at_gauss_pts_ptr, rho_n, rho_0));
278 fe_rhs_ptr->getOpPtrVector().push_back(
281 fe_rhs_ptr->getOpPtrVector().push_back(
282 new HookeElement::OpCalculateStrainAle(x_field, x_field,
284 fe_rhs_ptr->getOpPtrVector().push_back(
285 new HookeElement::OpCalculateStress<1>(x_field, x_field,
287 fe_rhs_ptr->getOpPtrVector().push_back(
288 new HookeElement::OpAleRhs_dx(x_field, x_field, data_at_pts));
289 fe_rhs_ptr->getOpPtrVector().push_back(
290 new HookeElement::OpCalculateEnergy(X_field, X_field,
292 fe_rhs_ptr->getOpPtrVector().push_back(
293 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
295 fe_rhs_ptr->getOpPtrVector().push_back(
296 new HookeElement::OpAleRhs_dX(X_field, X_field, data_at_pts));
301 CHKERR my_operators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
"x",
"X", ale,
304 CHKERR DMCreateGlobalVector(dm, &x);
305 CHKERR VecDuplicate(x, &f);
317 CHKERR MatDuplicate(
A, MAT_DO_NOT_COPY_VALUES, &fdA);
319 if (test_jacobian == PETSC_TRUE) {
320 char testing_options[] =
321 "-snes_test_jacobian -snes_test_jacobian_display "
322 "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 "
324 CHKERR PetscOptionsInsertString(NULL, testing_options);
326 char testing_options[] =
"-snes_no_convergence_test -snes_atol 0 "
327 "-snes_rtol 0 -snes_max_it 1 -pc_type none";
328 CHKERR PetscOptionsInsertString(NULL, testing_options);
332 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
337 CHKERR SNESSetFromOptions(snes);
339 CHKERR SNESSolve(snes, NULL, x);
341 if (test_jacobian == PETSC_FALSE) {
343 CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A0);
345 char testing_options_fd[] =
"-snes_fd";
346 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
350 CHKERR SNESSetFromOptions(snes);
352 CHKERR SNESSolve(snes, NULL, x);
353 CHKERR MatAXPY(
A, -1, fdA, SUBSET_NONZERO_PATTERN);
356 CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A);
357 PetscPrintf(PETSC_COMM_WORLD,
"Matrix norms %3.4e %3.4e\n", nrm_A,
361 const double tol = 1e-5;
364 "Difference between hand-calculated tangent matrix and finite "
365 "difference matrix is too big");
373 CHKERR SNESDestroy(&snes);
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
const std::string getDomainFEName() const
Get the Domain FE Name.
Interface for nonlinear (SNES) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
data for calculation heat conductivity and heat capacity elements
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &row_data)
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
OpGetDensityField(const std::string row_field, boost::shared_ptr< MatrixDouble > mat_coords_ptr, boost::shared_ptr< VectorDouble > density_at_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts_ptr)
boost::shared_ptr< MatrixDouble > matCoordsPtr
int operator()(int, int, int) const