19                                 {
   20 
   21  
   23 
   24  
   25  moab::Core mb_instance;              
   26  moab::Interface &moab = mb_instance; 
   27 
   28  try {
   29    
   32 
   33    PetscInt order_x = 2;
   34    PetscInt order_X = 2;
   35    PetscBool flg = PETSC_TRUE;
   36 
   37    PetscBool test_jacobian = PETSC_FALSE;
   39                               PETSC_NULLPTR);
   41                              &flg);
   43                              &flg);
   44 
   46 
   48 
   52                              3);
   54                                3);
   56 
   62 
   63    Range triangle_springs;
 
   65      if (it->getName().compare(0, 9, "SPRING_BC") == 0) {
   67                                                       triangle_springs, true);
   68      }
   69    }
   70 
   71    
   73                                           "MESH_NODE_POSITIONS");
   75        m_field, "SPATIAL_POSITION", "MESH_NODE_POSITIONS", triangle_springs);
   79 
   80    
   81    DM dm;
   83 
   84    PetscRandom rctx;
   85    PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
   86 
   87    auto set_coord = [&](
VectorAdaptor &&field_data, 
double *x, 
double *y,
 
   88                         double *z) {
   90      double value;
   92      PetscRandomGetValue(rctx, &value);
   93      field_data[0] = (*x) + (value - 0.5) * 
scale;
 
   94      PetscRandomGetValue(rctx, &value);
   95      field_data[1] = (*y) + (value - 0.5) * 
scale;
 
   96      PetscRandomGetValue(rctx, &value);
   97      field_data[2] = (*z) + (value - 0.5) * 
scale;
 
   99    };
  100 
  102                                                            "SPATIAL_POSITION");
  104        set_coord, "MESH_NODE_POSITIONS");
  105 
  106    PetscRandomDestroy(&rctx);
  107 
  108    boost::shared_ptr<NeumannForcesSurface> surfacePressure(
  110 
  111    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_rhs_ptr(
  112        surfacePressure, &(surfacePressure->getLoopFe()));
  113    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_lhs_ptr(
  114        surfacePressure, &(surfacePressure->getLoopFeLhs()));
  115    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_mat_rhs_ptr(
  116        surfacePressure, &(surfacePressure->getLoopFeMatRhs()));
  117    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE> fe_mat_lhs_ptr(
  118        surfacePressure, &(surfacePressure->getLoopFeMatLhs()));
  119 
  124 
  125    boost::shared_ptr<NeumannForcesSurface> surfaceForce(
  127 
  128    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
  129        fe_rhs_surface_force_ptr(surfaceForce, &(surfaceForce->getLoopFe()));
  130    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
  131        fe_lhs_surface_force_ptr(surfaceForce, &(surfaceForce->getLoopFeLhs()));
  132    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
  133        fe_mat_rhs_surface_force_ptr(surfaceForce,
  134                                     &(surfaceForce->getLoopFeMatRhs()));
  135    boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>
  136        fe_mat_lhs_surface_force_ptr(surfaceForce,
  137                                     &(surfaceForce->getLoopFeMatLhs()));
  138 
  140                          false, false);
  142                          false, false);
  144                          false, false);
  146                          false, false);
  147 
  149    CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes, 
false);
 
  150 
  151    nodes.pop_front();
  152    nodes.pop_back();
  153 
  154    boost::shared_ptr<NeumannForcesSurface::DataAtIntegrationPts> dataAtPts =
  155        boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
  156 
  157    dataAtPts->forcesOnlyOnEntitiesRow = nodes;
  158 
  160      if (
bit->getName().compare(0, 8, 
"PRESSURE") == 0) {
 
  161        CHKERR surfacePressure->addPressure(
"SPATIAL_POSITION", PETSC_NULLPTR,
 
  162                                            bit->getMeshsetId(), 
true, 
true);
 
  163        CHKERR surfacePressure->addPressureAle(
 
  164            "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
  166            true, true);
  167      }
  168    }
  169 
  170    const string block_set_force_name("FORCE");
  171    
  174      CHKERR surfaceForce->addForce(
"SPATIAL_POSITION", PETSC_NULLPTR,
 
  175                                    (
bit->getMeshsetId()), 
true, 
false);
 
  176      CHKERR surfaceForce->addForceAle(
 
  177          "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
  179          true, false, false);
  180    }
  181 
  183      if (it->getName().compare(0, block_set_force_name.length(),
  184                                block_set_force_name) == 0) {
  185        CHKERR surfaceForce->addForce(
"SPATIAL_POSITION", PETSC_NULLPTR,
 
  186                                      (it->getMeshsetId()), true, true);
  187        CHKERR surfaceForce->addForceAle(
 
  188            "SPATIAL_POSITION", "MESH_NODE_POSITIONS", dataAtPts,
  189            si->
getDomainFEName(), PETSC_NULLPTR, PETSC_NULLPTR, it->getMeshsetId(),
 
  190            true, true, false);
  191      }
  192    }
  193 
  194    
  195    
  196    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr(
  198    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr(
  201        m_field, fe_spring_lhs_ptr, fe_spring_rhs_ptr, "SPATIAL_POSITION",
  202        "MESH_NODE_POSITIONS");
  203 
  204    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ale_ptr_dx(
  206    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ale_ptr_dX(
  208 
  209    boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ale_ptr(
  211 
  212    boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
  213        data_at_spring_gp =
  214            boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(
  215                m_field);
  216 
  217    Range spring_ale_nodes;
 
  218    CHKERR moab.get_connectivity(triangle_springs, spring_ale_nodes, 
true);
 
  219 
  220    data_at_spring_gp->forcesOnlyOnEntitiesRow = spring_ale_nodes;
  221 
  223        m_field, fe_spring_lhs_ale_ptr_dx, fe_spring_lhs_ale_ptr_dX,
  224        fe_spring_rhs_ale_ptr, data_at_spring_gp, "SPATIAL_POSITION",
  226 
  228                                  PETSC_NULLPTR);
  230                                  PETSC_NULLPTR);
  231 
  233                                  PETSC_NULLPTR, PETSC_NULLPTR);
  235                                  PETSC_NULLPTR, PETSC_NULLPTR);
  237                                  PETSC_NULLPTR, PETSC_NULLPTR);
  238 
  240                                  nullptr, nullptr);
  242                                  nullptr, nullptr);
  243 
  245                                  nullptr, nullptr);
  247                                  nullptr, nullptr);
  248 
  249    
  251                                  fe_lhs_surface_force_ptr, nullptr, nullptr);
  253                                  fe_rhs_surface_force_ptr, nullptr, nullptr);
  254 
  256                                  fe_mat_lhs_surface_force_ptr, nullptr,
  257                                  nullptr);
  259                                  fe_mat_rhs_surface_force_ptr, nullptr,
  260                                  nullptr);
  261 
  263    CHKERR DMCreateGlobalVector(dm, &x);
 
  264    CHKERR VecDuplicate(x, &f);
 
  265    CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
 
  267 
  269    CHKERR DMCreateMatrix(dm, &A);
 
  270    CHKERR MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &fdA);
 
  271 
  272    if (test_jacobian == PETSC_TRUE) {
  273      char testing_options[] =
  274          "-snes_test_jacobian -snes_test_jacobian_display "
  275          "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 ";
  276      
  277      CHKERR PetscOptionsInsertString(NULL, testing_options);
 
  278    } else {
  279      char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
  280                               "-snes_rtol 0 "
  281                               "-snes_max_it 1 ";
  282      
  283      CHKERR PetscOptionsInsertString(NULL, testing_options);
 
  284    }
  285 
  286    SNES snes;
  287    CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
 
  292    CHKERR SNESSetFromOptions(snes);
 
  293 
  294    CHKERR SNESSolve(snes, NULL, x);
 
  295 
  296    if (test_jacobian == PETSC_FALSE) {
  297      double nrm_A0;
  298      CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
 
  299 
  300      char testing_options_fd[] = "-snes_fd";
  301      CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
 
  302 
  305      CHKERR SNESSetFromOptions(snes);
 
  306 
  307      CHKERR SNESSolve(snes, NULL, x);
 
  308      CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
 
  309 
  310      double nrm_A;
  311      CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
 
  312      PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
  313                  nrm_A / nrm_A0);
  314      nrm_A /= nrm_A0;
  315 
  316      const double tol = 1e-7;
 
  319                "Difference between hand-calculated tangent matrix and finite "
  320                "difference matrix is too big");
  321      }
  322    }
  323 
  328    CHKERR SNESDestroy(&snes);
 
  329 
  330    
  332  }
  334 
  335  
  337 
  338  return 0;
  339}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#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.
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
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
const FTensor::Tensor2< T, Dim, Dim > Vec
VectorShallowArrayAdaptor< double > VectorAdaptor
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 PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
virtual moab::Interface & get_moab()=0
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.
Simple interface for fast problem set-up.
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.
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
MoFEMErrorCode addBoundaryField(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 boundary.
std::vector< std::string > & getOtherFiniteElements()
Get the Other Finite Elements.
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 reference to pointer of interface.
Finite element and operators to apply force/pressures applied to surfaces.