diff --git a/CMakeLists.txt b/CMakeLists.txt index d09d3d66..27b5adcc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -179,6 +179,9 @@ set(scaleupROMObj_SOURCES include/dg_linear.hpp src/dg_linear.cpp + include/nlelast_integ.hpp + src/nlelast_integ.cpp + include/nonlinear_integ.hpp src/nonlinear_integ.cpp @@ -203,6 +206,9 @@ set(scaleupROMObj_SOURCES include/linelast_solver.hpp src/linelast_solver.cpp + include/nlelast_solver.hpp + src/nlelast_solver.cpp + include/stokes_solver.hpp src/stokes_solver.cpp diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp index 0ae4b356..d7ba2000 100644 --- a/include/linelast_solver.hpp +++ b/include/linelast_solver.hpp @@ -40,6 +40,7 @@ class LinElastSolver : public MultiBlockSolver Array as; // Lame constants for each subdomain, global boundary attribute ordering + double lambda, mu; Array lambda_c; Array mu_c; Array bdr_coeffs; @@ -53,7 +54,7 @@ class LinElastSolver : public MultiBlockSolver VectorFunctionCoefficient *init_x = NULL; public: - LinElastSolver(); + LinElastSolver(const double lambda_ = 1.0, const double mu_ = 1.0); virtual ~LinElastSolver(); diff --git a/include/nlelast_integ.hpp b/include/nlelast_integ.hpp new file mode 100644 index 00000000..aff0d423 --- /dev/null +++ b/include/nlelast_integ.hpp @@ -0,0 +1,200 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#ifndef SCALEUPROM_NLELAST_INTEG_HPP +#define SCALEUPROM_NLELAST_INTEG_HPP + +#include "mfem.hpp" +#include "hyperreduction_integ.hpp" + +namespace mfem +{ + + /// Abstract class for hyperelastic models that work with DG methods + class DGHyperelasticModel + { + protected: + ElementTransformation *Ttr; /**< Reference-element to target-element + transformation. */ + + public: + DGHyperelasticModel() : Ttr(NULL) {} + virtual ~DGHyperelasticModel() {} + + void SetTransformation(ElementTransformation &Ttr_) { Ttr = &Ttr_; } + virtual double EvalW(const DenseMatrix &Jpt) const = 0; + virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const = 0; + virtual void SetMatParam(ElementTransformation &Trans, const IntegrationPoint &ip)const = 0; + virtual void SetMatParam(FaceElementTransformations &Trans, const IntegrationPoint &ip)const = 0; + virtual void GetMatParam(Vector ¶ms) const = 0; + + virtual double EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) const = 0; + + virtual void EvalDmat(const int dim, const int dof, const DenseMatrix &DS, const DenseMatrix &J, DenseMatrix &Dmat)const = 0; + + }; + + class LinElastMaterialModel : public DGHyperelasticModel + { + protected: + Coefficient *c_mu, *c_lambda; + mutable double mu, lambda; + + public: + LinElastMaterialModel(double mu_, double lambda_) + : mu(mu_), lambda(lambda_) { c_mu = new ConstantCoefficient(mu), c_lambda = new ConstantCoefficient(lambda); } + + virtual void SetMatParam(ElementTransformation &Trans, const IntegrationPoint &ip) const; + virtual void SetMatParam(FaceElementTransformations &Trans, const IntegrationPoint &ip) const; + virtual void GetMatParam(Vector ¶ms) const; + virtual double EvalW(const DenseMatrix &J) const; + + virtual double EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) const; + virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const; + + virtual void EvalDmat(const int dim, const int dof, const DenseMatrix &DS, const DenseMatrix &J, DenseMatrix &Dmat) const; + + + }; + + class NeoHookeanHypModel : public DGHyperelasticModel + { + protected: + mutable double mu, K, g; + Coefficient *c_mu, *c_K, *c_g; + ElementTransformation *Ttr; + //DenseMatrix E, S; + mutable DenseMatrix Z; // dim x dim + mutable DenseMatrix G, C; // dof x dim + + public: + NeoHookeanHypModel(double mu_, double K_, double g_ = 1.0) + : mu(mu_), K(K_), g(g_) { c_mu = new ConstantCoefficient(mu), c_K = new ConstantCoefficient(K), c_g = new ConstantCoefficient(g); } + + virtual void SetMatParam(ElementTransformation &Trans, const IntegrationPoint &ip) const; + virtual void SetMatParam(FaceElementTransformations &Trans, const IntegrationPoint &ip) const; + virtual void GetMatParam(Vector ¶ms) const; + virtual double EvalW(const DenseMatrix &J) const; + + virtual double EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) const; + virtual void EvalP(const DenseMatrix &J, DenseMatrix &P) const; + + virtual void EvalDmat(const int dim, const int dof, const DenseMatrix &DS, const DenseMatrix &J, DenseMatrix &Dmat) const; + }; + + // DG boundary integrator for nonlinear elastic DG. + // For this is just DGElasticityIntegrator with a different name + class DGHyperelasticNLFIntegrator : virtual public HyperReductionIntegrator + { + + public: + DGHyperelasticNLFIntegrator(DGHyperelasticModel *m, double alpha_, double kappa_) + : HyperReductionIntegrator(false), model(m), lambda(NULL), mu(NULL), alpha(alpha_), kappa(kappa_) {} + + virtual void AssembleFaceVector(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Trans, + const Vector &elfun, Vector &elvect); + + virtual void AssembleFaceGrad(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Tr, + const Vector &elfun, DenseMatrix &elmat); + // values of all scalar basis functions for one component of u (which is a + protected: + DGHyperelasticModel *model; + Coefficient *lambda, *mu; + double alpha, kappa; // vector) at the integration point in the reference space + Vector elvect1, elvect2; + Vector elfun1, elfun2; + + DenseMatrix Jrt; + DenseMatrix PMatI1, PMatO1, DSh1, DS1, Jpt1, adjJ1, P1, Dmat1; + DenseMatrix PMatI2, PMatO2, DSh2, DS2, Jpt2, adjJ2, P2, Dmat2; + Vector shape1, tau1,wnor1; + Vector shape2, tau2,wnor2; + + Vector nor; // nor = |weight(J_face)| n + // 'jmat' corresponds to the term: kappa + DenseMatrix jmat; + + static void AssembleBlock( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, const Vector &row_shape, + const Vector &col_shape, const double jmatcoef, + const Vector &wnor, const DenseMatrix &Dmat, DenseMatrix &elmat, DenseMatrix &jmat); + + static void AssembleJmat( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, const Vector &row_shape, + const Vector &col_shape, const double jmatcoef, DenseMatrix &jmat); + }; + + // Domain integrator for nonlinear elastic DG. + class HyperelasticNLFIntegratorHR : virtual public HyperReductionIntegrator + { + + private: + DGHyperelasticModel *model; + // Jrt: the Jacobian of the target-to-reference-element transformation. + // Jpr: the Jacobian of the reference-to-physical-element transformation. + // Jpt: the Jacobian of the target-to-physical-element transformation. + // P: represents dW_d(Jtp) (dim x dim). + // DSh: gradients of reference shape functions (dof x dim). + // DS: gradients of the shape functions in the target (stress-free) + // configuration (dof x dim). + // PMatI: coordinates of the deformed configuration (dof x dim). + // PMatO: reshaped view into the local element contribution to the operator + // output - the result of AssembleElementVector() (dof x dim). + DenseMatrix DSh, DS, Jrt, Jpr, Jpt, P, PMatI, PMatO, Dmat; + + public: + HyperelasticNLFIntegratorHR(DGHyperelasticModel *m) : HyperReductionIntegrator(false), model(m) + { + } + + virtual void AssembleElementVector(const FiniteElement &el, + ElementTransformation &trans, + const Vector &elfun, + Vector &elvect); + + virtual void AssembleElementGrad(const FiniteElement &el, + ElementTransformation &trans, + const Vector &elfun, + DenseMatrix &elmat); + + virtual void AssembleH(const int dim, const int dof, const double w, + const DenseMatrix &J, DenseMatrix &elmat); + }; + + // RHS integrator for nonlinear elastic DG. + // For this is just DGElasticityDirichletLFIntegrator with a different name + class DGHyperelasticDirichletLFIntegrator : public LinearFormIntegrator // Should this be a nonlinear form later? + { + protected: + VectorCoefficient &uD; + DGHyperelasticModel *model; + Coefficient *lambda, *mu; + double alpha, kappa; + Vector shape, nor, u_dir; + + public: + DGHyperelasticDirichletLFIntegrator(VectorCoefficient &uD_, + DGHyperelasticModel *m, + double alpha_, double kappa_) + : uD(uD_), model(m), lambda(NULL), mu(NULL), alpha(alpha_), kappa(kappa_) {} + + virtual void AssembleRHSElementVect(const FiniteElement &el, + ElementTransformation &Tr, + Vector &elvect); + virtual void AssembleRHSElementVect(const FiniteElement &el, + FaceElementTransformations &Tr, + Vector &elvect); + + using LinearFormIntegrator::AssembleRHSElementVect; + }; + +} // namespace mfem + +#endif diff --git a/include/nlelast_solver.hpp b/include/nlelast_solver.hpp new file mode 100644 index 00000000..f8a5b914 --- /dev/null +++ b/include/nlelast_solver.hpp @@ -0,0 +1,154 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#ifndef SCALEUPROM_NLELAST_SOLVER_HPP +#define SCALEUPROM_NLELAST_SOLVER_HPP + +#include "multiblock_solver.hpp" +#include "interfaceinteg.hpp" +#include "mfem.hpp" +#include "nlelast_integ.hpp" + +// By convention we only use mfem namespace as default, not CAROM. +using namespace mfem; + + +// A proxy Operator used for FOM Newton Solver. +// Similar to SteadyNSOperator. +class NLElastOperator : public Operator +{ +protected: + bool direct_solve; + + mutable Vector x_u, y_u; + + Array u_offsets, vblock_offsets; + InterfaceForm *nl_itf = NULL; + Array hs; + BlockOperator *Hop = NULL; + + BlockMatrix *linearOp = NULL; + SparseMatrix *M = NULL, *B = NULL, *Bt = NULL; + + // Jacobian matrix objects + mutable BlockMatrix *system_jac = NULL; + mutable Array2D hs_mats; + mutable BlockMatrix *hs_jac = NULL; + mutable SparseMatrix *uu_mono = NULL; + mutable SparseMatrix *mono_jac = NULL; + mutable HypreParMatrix *jac_hypre = NULL; + + HYPRE_BigInt sys_glob_size; + mutable HYPRE_BigInt sys_row_starts[2]; +public: + NLElastOperator(const int height_, const int width_, Array &hs_, InterfaceForm *nl_itf_, + Array &u_offsets_, const bool direct_solve_=true); + + virtual ~NLElastOperator(); + + virtual void Mult(const Vector &x, Vector &y) const; + virtual Operator &GetGradient(const Vector &x) const; +}; + +class NLElastSolver : public MultiBlockSolver +{ + + friend class ParameterizedProblem; + friend class NLElastOperator; + +protected: + // interface integrator + InterfaceForm *a_itf = NULL; + // int skip_zeros = 1; + + // System matrix for Bilinear case. + Array2D mats; + // For nonlinear problem + // BlockOperator *globalMat; + BlockMatrix *globalMat = NULL; + SparseMatrix *globalMat_mono = NULL; + + // variables needed for direct solve + HYPRE_BigInt sys_glob_size; + HYPRE_BigInt sys_row_starts[2]; + HypreParMatrix *globalMat_hypre = NULL; + MUMPSSolver *mumps = NULL; + + // Temporary for nonlinear solve implementation + Array u_offsets, p_offsets, vblock_offsets; + BlockMatrix *systemOp = NULL; + Array hs; + InterfaceForm *nl_itf = NULL; + Solver *J_solver = NULL; + GMRESSolver *J_gmres = NULL; + NewtonSolver *newton_solver = NULL; + + // operators + Array bs; + Array as; + + // Lame constants for each subdomain, global boundary attribute ordering + Array lambda_c; + Array mu_c; + Array bdr_coeffs; + Array rhs_coeffs; + + // DG parameters specific to linear elasticity equation. + double alpha = -1.0; + double kappa = -1.0; + + DGHyperelasticModel* model; + + // Initial positions + VectorFunctionCoefficient *init_x = NULL; + +public: + NLElastSolver(DGHyperelasticModel* _model = NULL); + + virtual ~NLElastSolver(); + + static const std::vector GetVariableNames() + { + std::vector varnames(1); + varnames[0] = "solution"; + return varnames; + } + + virtual void InitVariables(); + + virtual void BuildOperators(); + virtual void BuildRHSOperators(); + virtual void BuildDomainOperators(); + + virtual void Assemble(); + virtual void AssembleRHS(); + virtual void AssembleOperator(); + // For bilinear case. + // system-specific. + virtual void AssembleInterfaceMatrices(); + + virtual bool Solve(); + + virtual void SetupBCVariables() override; + virtual void SetupIC(std::function F); + virtual void AddBCFunction(std::function F, const int battr = -1); + virtual void AddRHSFunction(std::function F); + virtual bool BCExistsOnBdr(const int &global_battr_idx); + virtual void SetupBCOperators(); + virtual void SetupRHSBCOperators(); + virtual void SetupDomainBCOperators(); + + // Component-wise assembly + virtual void BuildCompROMElement(Array &fes_comp); + virtual void BuildBdrROMElement(Array &fes_comp); + virtual void BuildInterfaceROMElement(Array &fes_comp); + + virtual void ProjectOperatorOnReducedBasis(); + + void SanityCheckOnCoeffs(); + + virtual void SetParameterizedProblem(ParameterizedProblem *problem); +}; + +#endif diff --git a/sketches/CMakeLists.txt b/sketches/CMakeLists.txt index e2ae78b7..c176ce9d 100644 --- a/sketches/CMakeLists.txt +++ b/sketches/CMakeLists.txt @@ -27,6 +27,10 @@ add_executable(precond precond.cpp $) add_executable(ns_mms ns_mms.cpp $) add_executable(ns_dg_mms ns_dg_mms.cpp $) add_executable(ns_rom ns_rom.cpp $) +add_executable(nlelast_mms nlelast_mms.cpp $) +add_executable(nlelast_cantilever nlelast_cantilever.cpp $) +add_executable(linelast_mms linelast_mms.cpp $) +add_executable(linelast_cantilever linelast_cantilever.cpp $) file(COPY inputs/gen_interface.yml DESTINATION ${CMAKE_BINARY_DIR}/sketches/inputs) file(COPY meshes/2x2.mesh DESTINATION ${CMAKE_BINARY_DIR}/sketches/meshes) diff --git a/sketches/linelast_mms.cpp b/sketches/linelast_mms.cpp new file mode 100644 index 00000000..ed72e83c --- /dev/null +++ b/sketches/linelast_mms.cpp @@ -0,0 +1,565 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "mfem.hpp" +#include +#include +#include +#include "linalg_utils.hpp" +#include "nonlinear_integ.hpp" +#include "nlelast_integ.hpp" +#include "dg_mixed_bilin.hpp" +#include "dg_bilinear.hpp" +#include "dg_linear.hpp" + +using namespace std; +using namespace mfem; + +static double K = 1.0; +static double mu = 1.0; + +// A proxy Operator used for FOM Newton Solver. +// Similar to SteadyNSOperator. +class SimpleNLElastOperator : public Operator +{ +protected: + NonlinearForm nlform; + + // Jacobian matrix objects + mutable SparseMatrix *J = NULL; + mutable HypreParMatrix *jac_hypre = NULL; + HYPRE_BigInt sys_glob_size; + mutable HYPRE_BigInt sys_row_starts[2]; + +public: + SimpleNLElastOperator(const int hw, NonlinearForm &nlform_); + + virtual ~SimpleNLElastOperator(); + + virtual void Mult(const Vector &x, Vector &y) const; + virtual Operator &GetGradient(const Vector &x) const; +}; + + +SimpleNLElastOperator::SimpleNLElastOperator(const int hw, NonlinearForm &nlform_) + : Operator(hw, hw),nlform(nlform_) +{ + // TODO: this needs to be changed for parallel implementation. + sys_glob_size = hw; + sys_row_starts[0] = 0; + sys_row_starts[1] = hw; +} + +SimpleNLElastOperator::~SimpleNLElastOperator() +{ + //delete J; + delete jac_hypre; + //delete nlform; +} + +void SimpleNLElastOperator::Mult(const Vector &x, Vector &y) const +{ + y = 0.0; + nlform.Mult(x, y); +} + +Operator& SimpleNLElastOperator::GetGradient(const Vector &x) const +{ + + J = dynamic_cast(&(nlform.GetGradient(x))); + jac_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, J); + return *jac_hypre; + //return nlform.GetGradient(x); +} + +// Define the analytical solution and forcing terms / boundary conditions +void SimpleExactRHSNeoHooke(const Vector &x, Vector &u); +void ExactSolutionNeoHooke(const Vector &x, Vector &u); +void SimpleExactSolutionNeoHooke(const Vector &x, Vector &u); +void ExactSolutionLinear(const Vector &x, Vector &u); +void ExactRHSLinear(const Vector &x, Vector &u); +void CheckGradient(NonlinearForm oper, FiniteElementSpace *fes); + + + +double error(Operator &M, Vector &x, Vector &b) +{ + assert(x.Size() == b.Size()); + + Vector res(x.Size()); + M.Mult(x, res); + res -= b; + + double tmp = 0.0; + for (int k = 0; k < x.Size(); k++) + tmp = max(tmp, abs(res(k))); + return tmp; +} + +double diff(Vector &a, Vector &b) +{ + assert(a.Size() == b.Size()); + + Vector res(a.Size()); + res = a; + res -= b; + + double tmp = 0.0; + for (int k = 0; k < a.Size(); k++) + tmp = max(tmp, abs(res(k))); + return tmp; +} + + + +int main(int argc, char *argv[]) +{ + int num_procs, rank; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + StopWatch chrono; + + // 1. Parse command-line options. + const char *mesh_file = "meshes/1x1.mesh"; + + int order = 1; + int refine = 0; + bool pa = false; + const char *device_config = "cpu"; + bool visualization = 1; + bool solve = false; + bool check_grad = false; + + OptionsParser args(argc, argv); + args.AddOption(&mesh_file, "-m", "--mesh", + "Mesh file to use."); + args.AddOption(&order, "-o", "--order", + "Finite element order (polynomial degree)."); + args.AddOption(&refine, "-r", "--refine", + "Number of refinements."); + args.AddOption(&solve, "-s", "--solve","-ns", "--nosolve", + "Solve the system."); + args.AddOption(&check_grad, "-cg", "--checkgrad","-ncg", "--nocheckgrad", + "Check gradients."); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(cout); + return 1; + } + args.PrintOptions(cout); + + // 3. Read the mesh from the given mesh file. We can handle triangular, + // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with + // the same code. + Mesh *mesh = new Mesh(mesh_file, 1, 1); + int dim = mesh->Dimension(); + + // 4. Refine the mesh to increase the resolution. In this example we do + // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the + // largest number that gives a final mesh with no more than 10,000 + // elements. + for (int l = 0; l < refine; l++) + { + mesh->UniformRefinement(); + } + + // 5. Define a finite element space on the mesh. Here we use the + FiniteElementCollection *dg_coll(new DG_FECollection(order, dim)); + FiniteElementCollection *h1_coll(new H1_FECollection(order, dim)); + + FiniteElementSpace *fes; + bool use_dg =true; + if (use_dg) + { + fes = new FiniteElementSpace(mesh, dg_coll, dim); + } + else + { + fes = new FiniteElementSpace(mesh, h1_coll, dim); + } + + // 7. Define the coefficients, analytical solution, and rhs of the PDE. + //VectorFunctionCoefficient exact_sol(dim, ExactSolutionNeoHooke); + //VectorFunctionCoefficient exact_RHS(dim, SimpleExactRHSNeoHooke); + VectorFunctionCoefficient exact_sol(dim, ExactSolutionLinear); + VectorFunctionCoefficient exact_RHS(dim, ExactRHSLinear); + + + // 8. Allocate memory (x, rhs) for the analytical solution and the right hand + // side. Define the GridFunction u,p for the finite element solution and + // linear forms fform and gform for the right hand side. The data + // allocated by x and rhs are passed as a reference to the grid functions + // (u,p) and the linear forms (fform, gform). + // MemoryType mt = device.GetMemoryType(); + int fomsize = fes->GetTrueVSize(); + + Vector x(fomsize), rhs(fomsize), rhs1(fomsize), rhs2(fomsize); + rhs = 0.0; + + // 12. Create the grid functions u and p. Compute the L2 error norms. + GridFunction u; + u.MakeRef(fes, x, 0); + + u = 0.0; + u.ProjectCoefficient(exact_sol); + + double kappa = (order+1)*(order+1); + //kappa = -1.0; + + //NeoHookeanHypModel model(mu, K); + LinElastMaterialModel model(mu, K); + + Array u_ess_attr(mesh->bdr_attributes.Max()); + // this array of integer essentially acts as the array of boolean: + // If value is 0, then it is not Dirichlet. + // If value is 1, then it is Dirichlet. + u_ess_attr = 1; + + Array u_ess_tdof; + fes->GetEssentialTrueDofs(u_ess_attr, u_ess_tdof); + + LinearForm *gform = new LinearForm(fes); + LinearForm *tempform = new LinearForm(fes); + gform->Update(fes, rhs1, 0); + gform->AddBdrFaceIntegrator(new DGHyperelasticDirichletLFIntegrator( + exact_sol, &model, 0.0, kappa), u_ess_attr); + gform->Assemble(); + gform->SyncAliasMemory(rhs1); + + cout<<"u_ess_tdof.Size() is: "<Update(fes, rhs2, 0); + tempform->AddDomainIntegrator(new VectorDomainLFIntegrator(exact_RHS)); + tempform->Assemble(); + gform->SyncAliasMemory(rhs2); + + + /* for (size_t i = 0; i < u_ess_tdof.Size(); i++) + { + rhs2[u_ess_tdof[i]] = 0.0; + } */ + + rhs += rhs1; + cout<<"rhs.Norml2() is: "<AddDomainIntegrator(new HyperelasticNLFIntegratorHR(&model)); + //nlform->AddBdrFaceIntegrator( new DGHyperelasticNLFIntegrator(&model, 0.0, kappa), u_ess_attr); + //nlform->AddBdrFaceIntegrator( new DGElasticityIntegrator(&model, 0.0, kappa), u_ess_attr); + //nlform->AddInteriorFaceIntegrator( new DGHyperelasticNLFIntegrator(&model, 0.0, kappa)); + //nlform->SetEssentialTrueDofs(u_ess_tdof); + + Vector lambda(mesh->attributes.Max()); + lambda = K; // Set lambda = 1 for all element attributes. + PWConstCoefficient lambda_c(lambda); + Vector _mu(mesh->attributes.Max()); + _mu = mu; // Set mu = 1 for all element attributes. + PWConstCoefficient mu_c(_mu); + + double alpha = -1.0; + + LinearForm b(fes); + b.AddBdrFaceIntegrator( + new DGElasticityDirichletLFIntegrator( + exact_sol, lambda_c, mu_c, alpha, kappa), u_ess_attr); + b.AddDomainIntegrator(new VectorDomainLFIntegrator(exact_RHS)); + b.Assemble(); + + BilinearForm a(fes); + a.AddDomainIntegrator(new ElasticityIntegrator(lambda_c, mu_c)); + a.AddBdrFaceIntegrator( + new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), u_ess_attr); + a.AddInteriorFaceIntegrator( + new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa)); + + // 10. Assemble the bilinear form and the corresponding linear system. + a.Assemble(); + + SparseMatrix A; + Vector B, X; + Array ess_tdof_list; + a.FormLinearSystem(ess_tdof_list, x, b, A, X, B); + + // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner. + //GSSmoother M(A); + //const double rtol = 1e-6; + //GMRES(A, M, B, X, 3, 5000, 1000, rtol*rtol, rtol*rtol); + + // 10. Recover the solution x as a grid function and save to file. The output + // can be viewed using GLVis as follows: "glvis -m mesh.mesh -g sol.gf" + + + SimpleNLElastOperator oper(fomsize, *nlform); + + if (check_grad) + { + CheckGradient(*nlform, fes); + /* SparseMatrix *jac = dynamic_cast(&(nlform->GetGradient(x))); + DenseMatrix J(*(jac->ToDenseMatrix())); + Vector ev(J.Size()); + cout<<"J.Det() is: "<GetNodes(x_ref); + //x_ref.ProjectCoefficient(exact_sol); + //_x = x_ref.GetTrueVector(); + + _y1 = 0.0; + oper.Mult(_x, _y1); //MFEM Neohookean + + double _y1_norm = _y1.Norml2(); + + + + /* for (size_t i = 0; i < _y1.Size(); i++) + { + cout<<"LHS(i) is: "<<_y1(i)< 0.001) + { + cout<<"LHS(i) is: "<<_y1[u_ess_tdof[i]]< 0.001) + { + cout<<"LHS(i) is: "<<_y1(i)<GetEssentialTrueDofs(ess_attr, ess_tdof); + + // 12. Create the grid functions u and p. Compute the L2 error norms. + GridFunction u(fes), us(fes); + for (int k = 0; k < u.Size(); k++) + { + u[k] = UniformRandom(); + us[k] = UniformRandom(); + } + + ConstantCoefficient one(1.0); + + Vector Nu(u.Size()); + oper.Mult(u, Nu); + double J0 = us * (Nu); + printf("J0: %.5E\n", J0); + + SparseMatrix *jac = dynamic_cast(&(oper.GetGradient(u))); + Vector grad(u.Size()); + jac->MultTranspose(us, grad); + double gg = grad * grad; + printf("gg: %.5E\n", gg); + + GridFunction u0(fes); + u0 = u; + + double error1 = 1.0e10; + printf("%10s\t%10s\t%10s\t%10s\t%10s\n", "amp", "J0", "J1", "dJdx", "error"); + for (int k = 0; k < 40; k++) + { + //double amp = pow(10.0, -0.25 * k); + double amp = pow(10.0, -5.0-0.25 * k); + double dx = amp; + if (gg > 1.0e-14) dx /= sqrt(gg); + + u.Set(1.0, u0); + u.Add(dx, grad); + + oper.Mult(u, Nu); + double J1 = us * (Nu); + double dJdx = (J1 - J0) / dx; + double error = abs((dJdx - gg)); + if (gg > 1.0e-14) error /= abs(gg); + + printf("%.5E\t%.5E\t%.5E\t%.5E\t%.5E\n", amp, J0, J1, dJdx, error); + + if (k > 4) + { + if (error > error1) + break; + else + error1 = error; + } + } + + return; +} diff --git a/sketches/nlelast_cantilever.cpp b/sketches/nlelast_cantilever.cpp new file mode 100644 index 00000000..ed5f7027 --- /dev/null +++ b/sketches/nlelast_cantilever.cpp @@ -0,0 +1,511 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "mfem.hpp" +#include +#include +#include +#include "linalg_utils.hpp" +#include "nonlinear_integ.hpp" +#include "nlelast_integ.hpp" +#include "dg_mixed_bilin.hpp" +#include "dg_bilinear.hpp" +#include "dg_linear.hpp" + +using namespace std; +using namespace mfem; + +static double K = 5.0; +static double mu = 5.0; +static const double pi = 4.0 * atan(1.0); +//static double f = -1.0; +static double f = 1.0; +//static double mu = 1.0; + +// A proxy Operator used for FOM Newton Solver. +// Similar to SteadyNSOperator. +class SimpleNLElastOperator : public Operator +{ +protected: + NonlinearForm nlform; + + // Jacobian matrix objects + mutable SparseMatrix *J = NULL; + mutable HypreParMatrix *jac_hypre = NULL; + HYPRE_BigInt sys_glob_size; + mutable HYPRE_BigInt sys_row_starts[2]; + +public: + SimpleNLElastOperator(const int hw, NonlinearForm &nlform_); + + virtual ~SimpleNLElastOperator(); + + virtual void Mult(const Vector &x, Vector &y) const; + virtual Operator &GetGradient(const Vector &x) const; +}; + + +SimpleNLElastOperator::SimpleNLElastOperator(const int hw, NonlinearForm &nlform_) + : Operator(hw, hw),nlform(nlform_) +{ + // TODO: this needs to be changed for parallel implementation. + sys_glob_size = hw; + sys_row_starts[0] = 0; + sys_row_starts[1] = hw; +} + +SimpleNLElastOperator::~SimpleNLElastOperator() +{ + //delete J; + delete jac_hypre; + //delete nlform; +} + +void SimpleNLElastOperator::Mult(const Vector &x, Vector &y) const +{ + y = 0.0; + nlform.Mult(x, y); +} + +Operator& SimpleNLElastOperator::GetGradient(const Vector &x) const +{ + + J = dynamic_cast(&(nlform.GetGradient(x))); + jac_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, J); + return *jac_hypre; + //return nlform.GetGradient(x); +} + +// Define the analytical solution and forcing terms / boundary conditions +void NullSolution(const Vector &x, Vector &u); +void NullDefSolution(const Vector &x, Vector &u); +void cantileverf(const Vector &x, Vector &u); +void CheckGradient(NonlinearForm oper, FiniteElementSpace *fes); + + +void CheckGradient(NonlinearForm oper, FiniteElementSpace *fes); + +double error(Operator &M, Vector &x, Vector &b) +{ + assert(x.Size() == b.Size()); + + Vector res(x.Size()); + M.Mult(x, res); + res -= b; + + double tmp = 0.0; + for (int k = 0; k < x.Size(); k++) + tmp = max(tmp, abs(res(k))); + return tmp; +} + +double diff(Vector &a, Vector &b) +{ + assert(a.Size() == b.Size()); + + Vector res(a.Size()); + res = a; + res -= b; + + double tmp = 0.0; + for (int k = 0; k < a.Size(); k++) + tmp = max(tmp, abs(res(k))); + return tmp; +} + + + +int main(int argc, char *argv[]) +{ + int num_procs, rank; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + StopWatch chrono; + + // 1. Parse command-line options. + const char *mesh_file = "meshes/1x1.mesh"; + + int order = 1; + int refine = 0; + bool pa = false; + const char *device_config = "cpu"; + bool visualization = 1; + bool solve = false; + bool check_grad = false; + bool nonlinear = true; + + OptionsParser args(argc, argv); + args.AddOption(&mesh_file, "-m", "--mesh", + "Mesh file to use."); + args.AddOption(&order, "-o", "--order", + "Finite element order (polynomial degree)."); + args.AddOption(&refine, "-r", "--refine", + "Number of refinements."); + args.AddOption(&solve, "-s", "--solve","-ns", "--nosolve", + "Solve the system."); + args.AddOption(&check_grad, "-cg", "--checkgrad","-ncg", "--nocheckgrad", + "Check gradients."); + args.AddOption(&nonlinear, "-nl", "--nonlinear","-lin", "--linear", + "toggle linear/nonlinear elasticity."); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(cout); + return 1; + } + args.PrintOptions(cout); + + // 3. Read the mesh from the given mesh file. We can handle triangular, + // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with + // the same code. + Mesh *mesh = new Mesh(mesh_file, 1, 1); + int dim = mesh->Dimension(); + + // 4. Refine the mesh to increase the resolution. In this example we do + // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the + // largest number that gives a final mesh with no more than 10,000 + // elements. + for (int l = 0; l < refine; l++) + { + mesh->UniformRefinement(); + } + + // 5. Define a finite element space on the mesh. Here we use the + FiniteElementCollection *dg_coll(new DG_FECollection(order, dim)); + FiniteElementCollection *h1_coll(new H1_FECollection(order, dim)); + + cout<<"dg_coll is: "<GetTrueVSize(); + + Vector x(fomsize), rhs(fomsize), rhs1(fomsize), rhs2(fomsize); + rhs = 0.0; + + // 12. Create the grid functions u and p. Compute the L2 error norms. + GridFunction u; + u.MakeRef(fes, x, 0); + + u = 0.0; + if (nonlinear) + { + u.ProjectCoefficient(*dir); + } + + + double kappa = (order+1)*(order+1); + + DGHyperelasticModel* model; + + if (nonlinear) + { +model = new NeoHookeanHypModel(mu, K); + } + else + { +model = new LinElastMaterialModel(mu, K); + } + + double alpha = -0.0; + + Array u_ess_attr(mesh->bdr_attributes.Max()); + // this array of integer essentially acts as the array of boolean: + // If value is 0, then it is not Dirichlet. + // If value is 1, then it is Dirichlet. + u_ess_attr = 0; + u_ess_attr[3] = 1; + + Array u_f_attr(mesh->bdr_attributes.Max()); + // this array of integer fentially acts as the array of boolean: + // If value is 0, then it is not Dirichlet. + // If value is 1, then it is Dirichlet. + u_f_attr = 0; + u_f_attr[1] = 1; + + Array u_ess_n_attr(mesh->bdr_attributes.Max()); + u_ess_attr = 1; + u_ess_attr[3] = 0; + u_ess_attr[1] = 0; + + Array u_ess_tdof; + fes->GetEssentialTrueDofs(u_ess_attr, u_ess_tdof); + + LinearForm *gform = new LinearForm(fes); + LinearForm *tempform = new LinearForm(fes); + gform->Update(fes, rhs1, 0); + gform->AddBdrFaceIntegrator(new DGHyperelasticDirichletLFIntegrator( + *dir, model, alpha, kappa), u_ess_attr); + gform->Assemble(); + gform->SyncAliasMemory(rhs1); + + cout<<"u_ess_tdof.Size() is: "<Update(fes, rhs2, 0); + tempform->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(*force), u_f_attr); + //tempform->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(*n_force), u_ess_n_attr); + //tempform->AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(*force), u_f_attr); + tempform->Assemble(); + gform->SyncAliasMemory(rhs2); + + + for (size_t i = 0; i < u_ess_tdof.Size(); i++) + { + rhs2[u_ess_tdof[i]] = 0.0; + } + + //rhs += rhs1; + cout<<"rhs.Norml2() is: "<AddDomainIntegrator(new HyperelasticNLFIntegratorHR(model)); + nlform->AddBdrFaceIntegrator(new DGHyperelasticNLFIntegrator(model, alpha, kappa), u_ess_attr); + //nlform->AddInteriorFaceIntegrator( new DGHyperelasticNLFIntegrator(model, 0.0, kappa)); + nlform->SetEssentialTrueDofs(u_ess_tdof); + + SimpleNLElastOperator oper(fomsize, *nlform); + + if (check_grad) + { + CheckGradient(*nlform, fes); + /* SparseMatrix *jac = dynamic_cast(&(nlform->GetGradient(x))); + DenseMatrix J(*(jac->ToDenseMatrix())); + Vector ev(J.Size()); + cout<<"J.Det() is: "<GetNodes(x_ref); + //x_ref.ProjectCoefficient(exact_sol); + //_x = x_ref.GetTrueVector(); + + _y1 = 0.0; + oper.Mult(_x, _y1); //MFEM Neohookean + + double _y1_norm = _y1.Norml2(); + + + + /* for (size_t i = 0; i < _y1.Size(); i++) + { + cout<<"LHS(i) is: "<<_y1(i)< 0.001) + { + cout<<"LHS(i) is: "<<_y1[u_ess_tdof[i]]< 0.001) + { + cout<<"LHS(i) is: "<<_y1(i)<GetEssentialTrueDofs(ess_attr, ess_tdof); + + // 12. Create the grid functions u and p. Compute the L2 error norms. + GridFunction u(fes), us(fes); + for (int k = 0; k < u.Size(); k++) + { + u[k] = UniformRandom(); + us[k] = UniformRandom(); + } + + ConstantCoefficient one(1.0); + + Vector Nu(u.Size()); + oper.Mult(u, Nu); + double J0 = us * (Nu); + printf("J0: %.5E\n", J0); + + SparseMatrix *jac = dynamic_cast(&(oper.GetGradient(u))); + Vector grad(u.Size()); + jac->MultTranspose(us, grad); + double gg = grad * grad; + printf("gg: %.5E\n", gg); + + GridFunction u0(fes); + u0 = u; + + double error1 = 1.0e10; + printf("%10s\t%10s\t%10s\t%10s\t%10s\n", "amp", "J0", "J1", "dJdx", "error"); + for (int k = 0; k < 40; k++) + { + //double amp = pow(10.0, -0.25 * k); + double amp = pow(10.0, -5.0-0.25 * k); + double dx = amp; + if (gg > 1.0e-14) dx /= sqrt(gg); + + u.Set(1.0, u0); + u.Add(dx, grad); + + oper.Mult(u, Nu); + double J1 = us * (Nu); + double dJdx = (J1 - J0) / dx; + double error = abs((dJdx - gg)); + if (gg > 1.0e-14) error /= abs(gg); + + printf("%.5E\t%.5E\t%.5E\t%.5E\t%.5E\n", amp, J0, J1, dJdx, error); + + if (k > 4) + { + if (error > error1) + break; + else + error1 = error; + } + } + + return; +} diff --git a/sketches/nlelast_mms.cpp b/sketches/nlelast_mms.cpp new file mode 100644 index 00000000..6b3ce598 --- /dev/null +++ b/sketches/nlelast_mms.cpp @@ -0,0 +1,542 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "mfem.hpp" +#include +#include +#include +#include "linalg_utils.hpp" +#include "nonlinear_integ.hpp" +#include "nlelast_integ.hpp" +#include "dg_mixed_bilin.hpp" +#include "dg_bilinear.hpp" +#include "dg_linear.hpp" + +using namespace std; +using namespace mfem; + +static double K = 1.0; +static double mu = 0.0; +static const double pi = 4.0 * atan(1.0); + +// A proxy Operator used for FOM Newton Solver. +// Similar to SteadyNSOperator. +class SimpleNLElastOperator : public Operator +{ +protected: + NonlinearForm nlform; + + // Jacobian matrix objects + mutable SparseMatrix *J = NULL; + mutable HypreParMatrix *jac_hypre = NULL; + HYPRE_BigInt sys_glob_size; + mutable HYPRE_BigInt sys_row_starts[2]; + +public: + SimpleNLElastOperator(const int hw, NonlinearForm &nlform_); + + virtual ~SimpleNLElastOperator(); + + virtual void Mult(const Vector &x, Vector &y) const; + virtual Operator &GetGradient(const Vector &x) const; +}; + + +SimpleNLElastOperator::SimpleNLElastOperator(const int hw, NonlinearForm &nlform_) + : Operator(hw, hw),nlform(nlform_) +{ + // TODO: this needs to be changed for parallel implementation. + sys_glob_size = hw; + sys_row_starts[0] = 0; + sys_row_starts[1] = hw; +} + +SimpleNLElastOperator::~SimpleNLElastOperator() +{ + //delete J; + delete jac_hypre; + //delete nlform; +} + +void SimpleNLElastOperator::Mult(const Vector &x, Vector &y) const +{ + y = 0.0; + nlform.Mult(x, y); +} + +Operator& SimpleNLElastOperator::GetGradient(const Vector &x) const +{ + + J = dynamic_cast(&(nlform.GetGradient(x))); + jac_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, J); + return *jac_hypre; + //return nlform.GetGradient(x); +} + +// Define the analytical solution and forcing terms / boundary conditions +void SimpleExactRHSNeoHooke(const Vector &x, Vector &u); +void ExactSolutionNeoHooke(const Vector &x, Vector &u); +void SimpleExactSolutionNeoHooke(const Vector &x, Vector &u); +void ExactSolutionLinear(const Vector &x, Vector &u); +void ExactRHSLinear(const Vector &x, Vector &u); +void CheckGradient(NonlinearForm oper, FiniteElementSpace *fes); + +double error(Operator &M, Vector &x, Vector &b) +{ + assert(x.Size() == b.Size()); + + Vector res(x.Size()); + M.Mult(x, res); + res -= b; + + double tmp = 0.0; + for (int k = 0; k < x.Size(); k++) + tmp = max(tmp, abs(res(k))); + return tmp; +} + +double diff(Vector &a, Vector &b) +{ + assert(a.Size() == b.Size()); + + Vector res(a.Size()); + res = a; + res -= b; + + double tmp = 0.0; + for (int k = 0; k < a.Size(); k++) + tmp = max(tmp, abs(res(k))); + return tmp; +} + + + +int main(int argc, char *argv[]) +{ + int num_procs, rank; + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &num_procs); + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + StopWatch chrono; + + // 1. Parse command-line options. + const char *mesh_file = "meshes/1x1.mesh"; + + int order = 1; + int refine = 0; + bool pa = false; + const char *device_config = "cpu"; + bool visualization = false; + bool solve = true; + bool check_grad = false; + bool nonlinear = true; + bool verbose = false; + bool weak_constraints = true; + + OptionsParser args(argc, argv); + args.AddOption(&mesh_file, "-m", "--mesh", + "Mesh file to use."); + args.AddOption(&order, "-o", "--order", + "Finite element order (polynomial degree)."); + args.AddOption(&refine, "-r", "--refine", + "Number of refinements."); + args.AddOption(&solve, "-s", "--solve","-ns", "--nosolve", + "Solve the system."); + args.AddOption(&check_grad, "-cg", "--checkgrad","-ncg", "--nocheckgrad", + "Check gradients."); + args.AddOption(&nonlinear, "-nl", "--nonlinear","-lin", "--linear", + "toggle linear/nonlinear elasticity."); + args.AddOption(&mu, "-mu", "--mu", + "Value of material parameter mu."); + args.AddOption(&K, "-K", "--K", + "Value of material parameter K."); + args.AddOption(&verbose, "-v", "--verbose","-nv", "--notverbose", + "Toggle verbose output."); + args.AddOption(&weak_constraints, "-wc", "--weakconstraints","-sc", "--strongconstraints", + "Toggle weak or strong constraint enforcement."); + args.AddOption(&visualization, "-vis", "--visualization","-nvis", "--novisualization", + "Toggle visualization."); + args.Parse(); + if (!args.Good()) + { + args.PrintUsage(cout); + return 1; + } + args.PrintOptions(cout); + + // 3. Read the mesh from the given mesh file. We can handle triangular, + // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with + // the same code. + Mesh *mesh = new Mesh(mesh_file, 1, 1); + int dim = mesh->Dimension(); + + // 4. Refine the mesh to increase the resolution. In this example we do + // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the + // largest number that gives a final mesh with no more than 10,000 + // elements. + for (int l = 0; l < refine; l++) + { + mesh->UniformRefinement(); + } + + // 5. Define a finite element space on the mesh. Here we use the + FiniteElementCollection *dg_coll(new DG_FECollection(order, dim)); + FiniteElementCollection *h1_coll(new H1_FECollection(order, dim)); + + FiniteElementSpace *fes; + bool use_dg = false; + if (use_dg) + { + fes = new FiniteElementSpace(mesh, dg_coll, dim); + } + else + { + fes = new FiniteElementSpace(mesh, h1_coll, dim); + } + + // 7. Define the coefficients, analytical solution, and rhs of the PDE. + VectorFunctionCoefficient* exact_sol; + VectorFunctionCoefficient* exact_RHS; + + if (nonlinear) + { + exact_sol = new VectorFunctionCoefficient(dim, ExactSolutionNeoHooke); + exact_RHS = new VectorFunctionCoefficient(dim, SimpleExactRHSNeoHooke); + } + else + { + mu = 1.0; + exact_sol = new VectorFunctionCoefficient(dim, ExactSolutionLinear); + exact_RHS = new VectorFunctionCoefficient(dim, ExactRHSLinear); + } + + // 8. Allocate memory (x, rhs) for the analytical solution and the right hand + // side. Define the GridFunction u,p for the finite element solution and + // linear forms fform and gform for the right hand side. The data + // allocated by x and rhs are passed as a reference to the grid functions + // (u,p) and the linear forms (fform, gform). + // MemoryType mt = device.GetMemoryType(); + int fomsize = fes->GetTrueVSize(); + + Vector x(fomsize), rhs(fomsize), rhs1(fomsize), rhs2(fomsize); + rhs = 0.0; + + // 12. Create the grid functions u and p. Compute the L2 error norms. + GridFunction u; + u.MakeRef(fes, x, 0); + + u = 0.0; + u.ProjectCoefficient(*exact_sol); + + double kappa = (order+1)*(order+1); + + DGHyperelasticModel* model; + if (nonlinear) + { + model = new NeoHookeanHypModel(mu, K); + } + else + { + model = new LinElastMaterialModel(mu, K); + } + + Array u_ess_attr(mesh->bdr_attributes.Max()); + Array u_ess_tdof; + u_ess_attr = 1; + + // Initialize the two RHS components + LinearForm *gform = new LinearForm(fes); + LinearForm *tempform = new LinearForm(fes); + + if (weak_constraints == false) + { + fes->GetEssentialTrueDofs(u_ess_attr, u_ess_tdof); + } + else + { + gform->Update(fes, rhs1, 0); + gform->AddBdrFaceIntegrator(new DGHyperelasticDirichletLFIntegrator( + *exact_sol, model, 0.0, kappa), u_ess_attr); + gform->Assemble(); + gform->SyncAliasMemory(rhs1); + rhs += rhs1; + if (verbose) + cout<<"rhs.Norml2() after added Dirichlet BCs is: "<Update(fes, rhs2, 0); + tempform->AddDomainIntegrator(new VectorDomainLFIntegrator(*exact_RHS)); + tempform->Assemble(); + tempform->SyncAliasMemory(rhs2); + + if (weak_constraints == false) + { + for (size_t i = 0; i < u_ess_tdof.Size(); i++) + { + rhs2[u_ess_tdof[i]] = 0.0; + } + kappa = 0.0; + } + + rhs += rhs2; + if (verbose) + cout<<"rhs.Norml2() after added forcing term is: "<AddDomainIntegrator(new HyperelasticNLFIntegratorHR(model)); + nlform->AddBdrFaceIntegrator( new DGHyperelasticNLFIntegrator(model, 0.0, kappa), u_ess_attr); + + if (weak_constraints == false) + { + nlform->SetEssentialTrueDofs(u_ess_tdof); + } + + SimpleNLElastOperator oper(fomsize, *nlform); + + if (check_grad) + { + CheckGradient(*nlform, fes); + } + + Vector _x(x); + + // Test applying nonlinear form + if (verbose) + { + Vector _y0(x); + Vector _y1(x); + + _y1 = 0.0; + oper.Mult(_x, _y1); + + double _y1_norm = _y1.Norml2(); + cout<<"--- RESIDUAL CHECK ---"< 0.1) + { + cout<<"High residual, res(i), is: "<GetEssentialTrueDofs(ess_attr, ess_tdof); + + // 12. Create the grid functions u and p. Compute the L2 error norms. + GridFunction u(fes), us(fes); + for (int k = 0; k < u.Size(); k++) + { + u[k] = UniformRandom(); + us[k] = UniformRandom(); + } + + ConstantCoefficient one(1.0); + + Vector Nu(u.Size()); + oper.Mult(u, Nu); + double J0 = us * (Nu); + printf("J0: %.5E\n", J0); + + SparseMatrix *jac = dynamic_cast(&(oper.GetGradient(u))); + Vector grad(u.Size()); + jac->MultTranspose(us, grad); + double gg = grad * grad; + printf("gg: %.5E\n", gg); + + GridFunction u0(fes); + u0 = u; + + double error1 = 1.0e10; + printf("%10s\t%10s\t%10s\t%10s\t%10s\n", "amp", "J0", "J1", "dJdx", "error"); + for (int k = 0; k < 40; k++) + { + //double amp = pow(10.0, -0.25 * k); + double amp = pow(10.0, -5.0-0.25 * k); + double dx = amp; + if (gg > 1.0e-14) dx /= sqrt(gg); + + u.Set(1.0, u0); + u.Add(dx, grad); + + oper.Mult(u, Nu); + double J1 = us * (Nu); + double dJdx = (J1 - J0) / dx; + double error = abs((dJdx - gg)); + if (gg > 1.0e-14) error /= abs(gg); + + printf("%.5E\t%.5E\t%.5E\t%.5E\t%.5E\n", amp, J0, J1, dJdx, error); + + if (k > 4) + { + if (error > error1) + break; + else + error1 = error; + } + } + + return; +} diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index a3e88e58..fc692b89 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -10,11 +10,15 @@ using namespace std; using namespace mfem; -LinElastSolver::LinElastSolver() +LinElastSolver::LinElastSolver(const double mu_, const double lambda_) : MultiBlockSolver() { alpha = config.GetOption("discretization/interface/alpha", -1.0); + alpha = 0.0; // TEMP kappa = config.GetOption("discretization/interface/kappa", (order + 1) * (order + 1)); + mu = mu_; + lambda = lambda_; + var_names = GetVariableNames(); num_var = var_names.size(); @@ -75,17 +79,7 @@ void LinElastSolver::SetupBCVariables() bdr_coeffs.SetSize(numBdr); bdr_coeffs = NULL; - lambda_c.SetSize(numSub); - lambda_c = NULL; - - mu_c.SetSize(numSub); - mu_c = NULL; - - for (size_t i = 0; i < numSub; i++) - { - lambda_c[i] = new ConstantCoefficient(1.0); - mu_c[i] = new ConstantCoefficient(1.0); - } + } void LinElastSolver::InitVariables() @@ -109,6 +103,18 @@ void LinElastSolver::InitVariables() var_offsets.PartialSum(); domain_offsets = var_offsets; + lambda_c.SetSize(numSub); + lambda_c = NULL; + + mu_c.SetSize(numSub); + mu_c = NULL; + + for (size_t i = 0; i < numSub; i++) + { + lambda_c[i] = new ConstantCoefficient(lambda); + mu_c[i] = new ConstantCoefficient(mu); + } + SetupBCVariables(); // Set up solution/rhs variables/ @@ -293,7 +299,15 @@ void LinElastSolver::AssembleOperator() globalMat_hypre = new HypreParMatrix(MPI_COMM_WORLD, sys_glob_size, sys_row_starts, globalMat_mono); mumps = new MUMPSSolver(MPI_COMM_SELF); + if (alpha != -1.0) + { + mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); + } + else + { mumps->SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE); + } + mumps->SetOperator(*globalMat_hypre); } } @@ -315,10 +329,14 @@ bool LinElastSolver::Solve() int print_level = config.GetOption("solver/print_level", 0); // TODO: need to change when the actual parallelization is implemented. - cout << "direct_solve is: " << direct_solve << endl; if (direct_solve) { assert(mumps); + if (alpha!=-1.0) + { + mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); + } + mumps->SetPrintLevel(print_level); mumps->Mult(*RHS, *U); } @@ -380,6 +398,11 @@ bool LinElastSolver::Solve() delete solver; } + PrintMatrix(*(globalMat->CreateMonolithic()->ToDenseMatrix()), "test_KLin.txt"); + PrintVector(*RHS, "RhsLin.txt"); + PrintVector(*U, "ULin2.txt"); + + return converged; } diff --git a/src/nlelast_integ.cpp b/src/nlelast_integ.cpp new file mode 100644 index 00000000..291b5a4e --- /dev/null +++ b/src/nlelast_integ.cpp @@ -0,0 +1,813 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "nlelast_integ.hpp" + +using namespace std; +namespace mfem +{ +double LinElastMaterialModel::EvalW(const DenseMatrix &J) const +{MFEM_ABORT("TODO")}; + +double LinElastMaterialModel::EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) const +{ +const double wL = w * c_lambda->Eval(Ttr, ip); +const double wM = w * c_mu->Eval(Ttr, ip); +return wL + 2.0*wM; +} + +void LinElastMaterialModel::SetMatParam(ElementTransformation &Trans, const IntegrationPoint &ip) const +{ + mu = c_mu->Eval(Trans, ip); + lambda = c_lambda->Eval(Trans, ip); +} + +void LinElastMaterialModel::SetMatParam(FaceElementTransformations &Trans, const IntegrationPoint &ip) const +{ + mu = c_mu->Eval(Trans, ip); + lambda = c_lambda->Eval(Trans, ip); +} + +void LinElastMaterialModel::GetMatParam(Vector ¶ms) const +{ + params.SetSize(2); + params[0] = lambda; + params[1] = mu; +} + +void LinElastMaterialModel::EvalP(const DenseMatrix &J, DenseMatrix &P) const +{ + // stress = 2*M*e(u) + L*tr(e(u))*I, where + // e(u) = (1/2)*(grad(u) + grad(u)^T) + const int dim = J.Size(); + double L = lambda; + const double mu2 = 2.0*mu; + if (dim == 2) + { + L *= (J(0,0) + J(1,1)); + P(0,0) = mu2*J(0,0) + L; + P(1,1) = mu2*J(1,1) + L; + P(1,0) = mu*(J(0,1) + J(1,0)); + P(0,1) = mu*(J(0,1) + J(1,0)); + } + else if (dim == 3) + { + L *= (J(0,0) + J(1,1) + J(2,2)); + P(0,0) = mu2*J(0,0) + L; + P(1,1) = mu2*J(1,1) + L; + P(2,2) = mu2*J(2,2) + L; + P(0,1) = mu*(J(0,1) + J(1,0)); + P(1,0) = mu*(J(0,1) + J(1,0)); + P(0,2) = mu*(J(0,2) + J(2,0)); + P(2,0) = mu*(J(0,2) + J(2,0)); + P(1,2) = mu*(J(1,2) + J(2,1)); + P(2,1) = mu*(J(1,2) + J(2,1)); + } +} + +void LinElastMaterialModel::EvalDmat(const int dim, const int dof, const DenseMatrix &DS, const DenseMatrix &J, DenseMatrix &Dmat) const +{ +for (size_t i = 0; i < dim; i++) + { + for (size_t j = 0; j < dim; j++) // Looping over each entry in residual + { + const int S_ij = j * dim + i; + + for (size_t m = 0; m < dof; m++) + for (size_t n = 0; n < dim; n++) // Looping over derivatives with respect to U + { + const int U_mn = n * dof + m; + Dmat(S_ij, U_mn) = ((i == j) ? lambda * DS(m,n) : 0.0); + Dmat(S_ij, U_mn) += ((n == i) ? mu * DS(m,j) : 0.0); + Dmat(S_ij, U_mn) += ((n == j) ? mu * DS(m,i) : 0.0); + } + } + } +} + +double NeoHookeanHypModel::EvalW(const DenseMatrix &J) const +{ int dim = J.Width(); + + double dJ = J.Det(); + double sJ = dJ/g; + double bI1 = pow(dJ, -2.0/dim)*(J*J); // \bar{I}_1 + + return 0.5*(mu*(bI1 - dim) + K*(sJ - 1.0)*(sJ - 1.0)); + } + +double NeoHookeanHypModel::EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) const +{ +const double wL = w * c_K->Eval(Ttr, ip); +const double wM = w * c_mu->Eval(Ttr, ip); +return wL + 2.0*wM; +} + +void NeoHookeanHypModel::SetMatParam(ElementTransformation &Trans, const IntegrationPoint &ip) const +{ + mu = c_mu->Eval(Trans, ip); + K = c_K->Eval(Trans, ip); +} + +void NeoHookeanHypModel::SetMatParam(FaceElementTransformations &Trans, const IntegrationPoint &ip) const +{ + mu = c_mu->Eval(Trans, ip); + K = c_K->Eval(Trans, ip); +} + +void NeoHookeanHypModel::GetMatParam(Vector ¶ms) const +{ + params.SetSize(2); + params[0] = K; + params[1] = mu; +} + +void NeoHookeanHypModel::EvalP(const DenseMatrix &J, DenseMatrix &P) const +{ + int dim = J.Width(); + + Z.SetSize(dim); + CalcAdjugateTranspose(J, Z); + + double dJ = J.Det(); + double a = mu*pow(dJ, -2.0/dim); + double b = K*(dJ/g - 1.0)/g - a*(J*J)/(dim*dJ); + + P = 0.0; + P.Add(a, J); + P.Add(b, Z); +} + +void NeoHookeanHypModel::EvalDmat(const int dim, const int dof, const DenseMatrix &DS, const DenseMatrix &J, DenseMatrix &Dmat) const +{ + double dJ = J.Det(); + double sJ = dJ/g; + double a = mu*pow(dJ, -2.0/dim); + double bc = a*(J*J)/dim; + double b = bc - K*sJ*(sJ - 1.0); + double c = 2.0*bc/dim + K*sJ*(2.0*sJ - 1.0); + + Z.SetSize(dim); + G.SetSize(dof, dim); + C.SetSize(dof, dim); + + CalcAdjugateTranspose(J, Z); + Z *= (1.0/dJ); // Z = J^{-t} + + MultABt(DS, J, C); // C = DS J^t + MultABt(DS, Z, G); // G = DS J^{-1} + + const double a2 = a * (-2.0/dim); + + for (size_t i = 0; i < dim; i++) + for (size_t j = 0; j < dim; j++) // Looping over each entry in residual + { + const int S_ij = j * dim + i; + + for (size_t m = 0; m < dof; m++) + for (size_t n = 0; n < dim; n++) // Looping over derivatives with respect to U + { + const int U_mn = n * dof + m; + + const double s1 = (i==n) ? a * DS(m,j) : 0.0; + const double s2 = a2 * (J(i,j)*G(m,n) + Z(i,j)*C(m,n)) + + b*Z(n,j)*G(m,i) + c*Z(i,j)*G(m,n); + + Dmat(S_ij, U_mn) = s1 + s2; + } + } + +} + +void _PrintMatrix(const DenseMatrix &mat, + const std::string &filename) +{ + FILE *fp = fopen(filename.c_str(), "w"); + + for (int i = 0; i < mat.NumRows(); i++) + { + for (int j = 0; j < mat.NumCols(); j++) + fprintf(fp, "%.15E\t", mat(i,j)); + fprintf(fp, "\n"); + } + + fclose(fp); + return; +} + +void _PrintVector(const Vector &vec, + const std::string &filename) +{ + FILE *fp = fopen(filename.c_str(), "w"); + + for (int i = 0; i < vec.Size(); i++) + fprintf(fp, "%.15E\n", vec(i)); + + fclose(fp); + return; +} + +void DGHyperelasticNLFIntegrator::AssembleJmat( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, const Vector &row_shape, + const Vector &col_shape, const double jmatcoef,DenseMatrix &jmat){ + for (int d = 0; d < dim; ++d) + { + const int jo = col_offset + d*col_ndofs; + const int io = row_offset + d*row_ndofs; + for (int jdof = 0, j = jo; jdof < col_ndofs; ++jdof, ++j) + { + const double sj = jmatcoef * col_shape(jdof); + for (int i = max(io,j), idof = i - io; idof < row_ndofs; ++idof, ++i) + { + jmat(i, j) += row_shape(idof) * sj; + } + } + } +}; + + // Boundary integrator +void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Trans, + const Vector &elfun, Vector &elvect) +{ + const int dim = el1.GetDim(); + const int ndofs1 = el1.GetDof(); + const int ndofs2 = (Trans.Elem2No >= 0) ? el2.GetDof() : 0; + + const int nvdofs = dim*(ndofs1 + ndofs2); + + Vector elfun_copy(elfun); // FIXME: How to avoid this? + + nor.SetSize(dim); + Jrt.SetSize(dim); + elvect.SetSize(nvdofs); + elvect = 0.0; + + const bool kappa_is_nonzero = (kappa != 0.0); + if (kappa_is_nonzero) + { + jmat.SetSize(nvdofs); + jmat = 0.; + } + + model->SetTransformation(Trans); + + shape1.SetSize(ndofs1); + elfun1.MakeRef(elfun_copy,0,ndofs1*dim); + elvect1.MakeRef(elvect,0,ndofs1*dim); + PMatI1.UseExternalData(elfun1.GetData(), ndofs1, dim); + DSh1.SetSize(ndofs1, dim); + DS1.SetSize(ndofs1, dim); + Jpt1.SetSize(dim); + adjJ1.SetSize(dim); + P1.SetSize(dim); + tau1.SetSize(dim); + + if (ndofs2) + { + shape2.SetSize(ndofs2); + elfun2.MakeRef(elfun_copy,ndofs1*dim,ndofs2*dim); + elvect2.MakeRef(elvect,ndofs1*dim,ndofs2*dim); + PMatI2.UseExternalData(elfun2.GetData(), ndofs2, dim); + DSh2.SetSize(ndofs2, dim); + DS2.SetSize(ndofs2, dim); + adjJ2.SetSize(dim); + Jpt2.SetSize(dim); + P2.SetSize(dim); + tau2.SetSize(dim); + } + + const IntegrationRule *ir = IntRule; + if (ir == NULL) + { + int intorder = 2*el1.GetOrder(); + ir = &IntRules.Get(Trans.GetGeometryType(), intorder); + } + + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + + // Set the integration point in the face and the neighboring element + Trans.SetAllIntPoints(&ip); + + // Access the neighboring element's integration point + const IntegrationPoint &eip1 = Trans.GetElement1IntPoint(); + const IntegrationPoint &eip2 = Trans.GetElement2IntPoint(); + + if (dim == 1) + { + nor(0) = 2*eip1.x - 1.0; + } + else + { + CalcOrtho(Trans.Jacobian(), nor); + } + + CalcInverse(Trans.Jacobian(), Jrt); + + double w = ip.weight; + double wLM = 0.0; + if (ndofs2) + { + w /= 2.0; + + el2.CalcShape(eip2, shape2); + el2.CalcDShape(eip2, DSh2); + CalcAdjugate(Trans.Elem2->Jacobian(), adjJ2); + Mult(DSh2, adjJ2, DS2); + MultAtB(PMatI2, DS2, Jpt2); + model->SetMatParam(Trans, eip2); + + model->EvalP(Jpt2, P2); + + double w2 = w / Trans.Elem2->Weight(); + wLM = model->EvalDGWeight(w2, *Trans.Elem2, eip2); + P2 *= w2; + P2.Mult(nor, tau2); + } + + el1.CalcShape(eip1, shape1); + el1.CalcDShape(eip1, DSh1); + CalcAdjugate(Trans.Elem1->Jacobian(), adjJ1); + Mult(DSh1, adjJ1, DS1); + MultAtB(PMatI1, DS1, Jpt1); + model->SetMatParam(Trans, eip1); + + model->EvalP(Jpt1, P1); + + double w1 = w / Trans.Elem1->Weight(); + wLM += model->EvalDGWeight(w1, *Trans.Elem1, eip1); + P1 *= w1; + + P1.Mult(nor, tau1); + + const double jmatcoef = kappa * (nor*nor) * wLM; + + // (1,1) block + for (int im = 0, i = 0; im < dim; ++im) + { + for (int idof = 0; idof < ndofs1; ++idof, ++i) + { + elvect(i) += shape1(idof) * tau1(im); + } + } + if (ndofs2 != 0) { + shape2.Neg(); + } + + if (kappa_is_nonzero) + { + jmat = 0.; + AssembleJmat( + dim, ndofs1, ndofs1, 0, 0, shape1, + shape1, jmatcoef, jmat); + if (ndofs2 != 0) { + AssembleJmat( + dim, ndofs1, ndofs2, 0, ndofs1*dim, shape1, + shape2, jmatcoef, jmat); + AssembleJmat( + dim, ndofs2, ndofs1, ndofs1*dim, 0, shape2, + shape1, jmatcoef, jmat); + AssembleJmat( + dim, ndofs2, ndofs2, ndofs1*dim, ndofs1*dim, shape2, + shape2, jmatcoef, jmat); + } + + //Flatten jmat + for (int i = 0; i < nvdofs; ++i) + { + for (int j = 0; j < i; ++j) + { + jmat(j,i) = jmat(i,j); + } + } + // Apply jmat + for (size_t i = 0; i < nvdofs; i++) + { + for (size_t j = 0; j < nvdofs; j++) + { + elvect(i) -= jmat(i,j) * (elfun(j)); + } + } + + } + + if (ndofs2 == 0) {continue;} + + // (1,2) block + for (int im = 0, i = 0; im < dim; ++im) + { + for (int idof = 0; idof < ndofs1; ++idof, ++i) + { + elvect(i) += shape1(idof) * tau2(im); + } + } + + // (2,1) block + for (int im = 0, i = ndofs1*dim; im < dim; ++im) + { + for (int idof = 0; idof < ndofs2; ++idof, ++i) + { + elvect(i) += shape2(idof) * tau1(im); + } + } + + // (2,2) block + for (int im = 0, i = ndofs1*dim; im < dim; ++im) + { + for (int idof = 0; idof < ndofs2; ++idof, ++i) + { + elvect(i) += shape2(idof) * tau2(im); + } + } + + } + + elvect *= -1.0; +} + +void DGHyperelasticNLFIntegrator::AssembleFaceGrad(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Tr, + const Vector &elfun, DenseMatrix &elmat){ + +const int dim = el1.GetDim(); + const int ndofs1 = el1.GetDof(); + const int ndofs2 = (Tr.Elem2No >= 0) ? el2.GetDof() : 0; + + const int nvdofs = dim*(ndofs1 + ndofs2); + + Vector elfun_copy(elfun); // FIXME: How to avoid this? + nor.SetSize(dim); + Jrt.SetSize(dim); + + elmat.SetSize(nvdofs); + + elmat = 0.0; + model->SetTransformation(Tr); + + const bool kappa_is_nonzero = (kappa != 0.0); + if (kappa_is_nonzero) + { + jmat.SetSize(nvdofs); + jmat = 0.; + } + + shape1.SetSize(ndofs1); + elfun1.MakeRef(elfun_copy,0,ndofs1*dim); + PMatI1.UseExternalData(elfun1.GetData(), ndofs1, dim); + DSh1.SetSize(ndofs1, dim); + DS1.SetSize(ndofs1, dim); + Jpt1.SetSize(dim); + P1.SetSize(dim); + tau1.SetSize(dim); + wnor1.SetSize(dim); + adjJ1.SetSize(dim); + Dmat1.SetSize(dim * dim, dim*ndofs1); + Dmat1 = 0.0; + + if (ndofs2) + { + shape2.SetSize(ndofs2); + elfun2.MakeRef(elfun_copy,ndofs1*dim,ndofs2*dim); + PMatI2.UseExternalData(elfun2.GetData(), ndofs2, dim); + DSh2.SetSize(ndofs2, dim); + DS2.SetSize(ndofs2, dim); + Jpt2.SetSize(dim); + P2.SetSize(dim); + tau2.SetSize(dim); + wnor2.SetSize(dim); + adjJ2.SetSize(dim); + Dmat2.SetSize(dim * dim, dim*ndofs2); + Dmat2 = 0.0; + } + + const IntegrationRule *ir = IntRule; + if (ir == NULL) + { + int intorder = 2*el1.GetOrder(); + ir = &IntRules.Get(Tr.GetGeometryType(), intorder); + } + + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + + // Set the integration point in the face and the neighboring element + Tr.SetAllIntPoints(&ip); + + // Access the neighboring element's integration point + const IntegrationPoint &eip1 = Tr.GetElement1IntPoint(); + const IntegrationPoint &eip2 = Tr.GetElement2IntPoint(); + + if (dim == 1) + { + nor(0) = 2*eip1.x - 1.0; + } + else + { + CalcOrtho(Tr.Jacobian(), nor); + } + + CalcInverse(Tr.Jacobian(), Jrt); + + double w = ip.weight; + double wLM = 0.0; + if (ndofs2) + { + w /= 2.0; + el2.CalcShape(eip2, shape2); + el2.CalcDShape(eip2, DSh2); + + CalcAdjugate(Tr.Elem2->Jacobian(), adjJ2); + Mult(DSh2, adjJ2, DS2); + + model->SetMatParam(Tr,eip2); + + MultAtB(PMatI2, DS2, Jpt2); + + model->EvalDmat(dim, ndofs2, DS2, Jpt2, Dmat2); + double w2 = w / Tr.Elem2->Weight(); + wLM = model->EvalDGWeight(w2, *Tr.Elem2, eip2); + wnor2.Set(w2,nor); + } + + el1.CalcShape(eip1, shape1); + el1.CalcDShape(eip1, DSh1); + + double w1 = w / Tr.Elem1->Weight(); + wLM += model->EvalDGWeight(w1, *Tr.Elem1, eip1); + + CalcAdjugate(Tr.Elem1->Jacobian(), adjJ1); + Mult(DSh1, adjJ1, DS1); + + model->SetMatParam(Tr,eip1); + MultAtB(PMatI1, DS1, Jpt1); + model->EvalDmat(dim, ndofs1, DS1, Jpt1, Dmat1); + + const double jmatcoef = kappa * (nor*nor) * wLM; + + // (1,1) block + wnor1.Set(w1,nor); + AssembleBlock(dim, ndofs1, ndofs1, 0, 0, shape1, shape1, jmatcoef, wnor1,Dmat1, elmat,jmat); + + if (ndofs2 == 0) {continue;} + shape2.Neg(); + + // (1,2) block + AssembleBlock(dim, ndofs1, ndofs2, 0, dim*ndofs1, shape1, shape2,jmatcoef,wnor2, Dmat2, elmat,jmat); + + // (2,1) block + AssembleBlock(dim, ndofs2, ndofs1, dim*ndofs1, 0, shape2, shape1,jmatcoef,wnor1, Dmat1, elmat,jmat); + + // (2,2) block + AssembleBlock(dim, ndofs2, ndofs2, dim*ndofs1, dim*ndofs1, shape2, shape2,jmatcoef,wnor2, Dmat2, elmat,jmat); + } + + // elmat := -elmat + jmat + elmat *= -1.0; + if (kappa_is_nonzero) + { + for (int i = 0; i < nvdofs; ++i) + { + for (int j = 0; j < i; ++j) + { + double mij = jmat(i,j); + elmat(i,j) += mij; + elmat(j,i) += mij; + } + elmat(i,i) += jmat(i,i); + } + } + + }; + +void DGHyperelasticNLFIntegrator::AssembleBlock( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, const Vector &row_shape, + const Vector &col_shape, const double jmatcoef, + const Vector &wnor, const DenseMatrix &Dmat, DenseMatrix &elmat, DenseMatrix &jmat){ +for (int n = 0, jj = col_offset; n < dim; ++n) + { + for (int m = 0; m < col_ndofs; ++m, ++jj) + { + const int mn = n * col_ndofs + m; + for (int j = 0, ii = row_offset; j < dim; ++j) + { + for (int i = 0; i < row_ndofs; ++i, ++ii) + { + double temp = 0.0; + for (size_t k = 0; k < dim; k++) + { + const int S_jk = k * dim + j; + temp += Dmat(S_jk, mn) * wnor(k); + } + elmat(ii, jj) += row_shape(i) * temp; + } + } + } + } + + if (jmatcoef == 0.0) { return; } + + for (int d = 0; d < dim; ++d) + { + const int jo = col_offset + d*col_ndofs; + const int io = row_offset + d*row_ndofs; + for (int jdof = 0, j = jo; jdof < col_ndofs; ++jdof, ++j) + { + const double sj = jmatcoef * col_shape(jdof); + for (int i = max(io,j), idof = i - io; idof < row_ndofs; ++idof, ++i) + { + jmat(i, j) += row_shape(idof) * sj; + } + } + } +}; + +// Domain integrator +void HyperelasticNLFIntegratorHR::AssembleElementVector(const FiniteElement &el, + ElementTransformation &Ttr, + const Vector &elfun, + Vector &elvect){ + int dof = el.GetDof(), dim = el.GetDim(); + + DSh.SetSize(dof, dim); + DS.SetSize(dof, dim); + Jrt.SetSize(dim); + Jpt.SetSize(dim); + P.SetSize(dim); + PMatI.UseExternalData(elfun.GetData(), dof, dim); + elvect.SetSize(dof*dim); + PMatO.UseExternalData(elvect.GetData(), dof, dim); + + const IntegrationRule *ir = IntRule; + if (!ir) + { + ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <--- + } + + elvect = 0.0; + model->SetTransformation(Ttr); + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + Ttr.SetIntPoint(&ip); + CalcInverse(Ttr.Jacobian(), Jrt); + + el.CalcDShape(ip, DSh); + Mult(DSh, Jrt, DS); + MultAtB(PMatI, DS, Jpt); + + model->EvalP(Jpt, P); + + P *= ip.weight * Ttr.Weight(); + AddMultABt(DS, P, PMatO); + } + } + +void HyperelasticNLFIntegratorHR::AssembleElementGrad(const FiniteElement &el, + ElementTransformation &Ttr, + const Vector &elfun, + DenseMatrix &elmat){ + int dof = el.GetDof(), dim = el.GetDim(); + Dmat.SetSize(dim * dim, dim * dof); + + DSh.SetSize(dof, dim); + DS.SetSize(dof, dim); + Jrt.SetSize(dim); + Jpt.SetSize(dim); + PMatI.UseExternalData(elfun.GetData(), dof, dim); + elmat.SetSize(dof*dim); + + const IntegrationRule *ir = IntRule; + if (!ir) + { + ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <--- + } + + elmat = 0.0; + model->SetTransformation(Ttr); + + for (int i = 0; i < ir->GetNPoints(); i++) + { + const IntegrationPoint &ip = ir->IntPoint(i); + Ttr.SetIntPoint(&ip); + CalcInverse(Ttr.Jacobian(), Jrt); + + el.CalcDShape(ip, DSh); + Mult(DSh, Jrt, DS); + MultAtB(PMatI, DS, Jpt); + + model->SetMatParam(Ttr,ip); + + model->EvalDmat(dim, dof, DS, Jpt, Dmat); + AssembleH(dim, dof, ip.weight * Ttr.Weight(), DS, elmat); + } + }; + + +void HyperelasticNLFIntegratorHR::AssembleH(const int dim, const int dof, const double w, const DenseMatrix &J, DenseMatrix &elmat) + { + //Assemble elmat + for (size_t i = 0; i < dof; i++) + { + for (size_t j = 0; j < dim; j++) // Looping over each entry in residual + { + const int ij = j * dof + i; + + for (size_t m = 0; m < dof; m++) + for (size_t n = 0; n < dim; n++) // Looping over derivatives with respect to U + { + const int mn = n * dof + m; + double temp = 0.0; + for (size_t k = 0; k < dim; k++) + { + const int S_jk = k * dim + j; + temp += Dmat(S_jk, mn) * w * J(i,k); + } + elmat(ij, mn) += temp; + } + } + } + + }; + + +//RHS +void DGHyperelasticDirichletLFIntegrator::AssembleRHSElementVect(const FiniteElement &el, + ElementTransformation &Tr, + Vector &elvect){ + mfem_error("DGElasticityDirichletLFIntegrator::AssembleRHSElementVect"); + }; + +void DGHyperelasticDirichletLFIntegrator::AssembleRHSElementVect(const FiniteElement &el, + FaceElementTransformations &Tr, + Vector &elvect){ +MFEM_ASSERT(Tr.Elem2No < 0, "interior boundary is not supported"); + + const int dim = el.GetDim(); + const int ndofs = el.GetDof(); + const int nvdofs = dim*ndofs; + + elvect.SetSize(nvdofs); + elvect = 0.0; + + shape.SetSize(ndofs); + nor.SetSize(dim); + u_dir.SetSize(dim); + + const IntegrationRule *ir = IntRule; + if (ir == NULL) + { + const int order = 2*el.GetOrder(); // <----- + ir = &IntRules.Get(Tr.GetGeometryType(), order); + } + + for (int pi = 0; pi < ir->GetNPoints(); ++pi) + { + const IntegrationPoint &ip = ir->IntPoint(pi); + + // Set the integration point in the face and the neighboring element + Tr.SetAllIntPoints(&ip); + + // Access the neighboring element's integration point + const IntegrationPoint &eip = Tr.GetElement1IntPoint(); + + // Evaluate the Dirichlet b.c. using the face transformation. + uD.Eval(u_dir, Tr, ip); + el.CalcShape(eip, shape); + + if (dim == 1) + { + nor(0) = 2*eip.x - 1.0; + } + else + { + CalcOrtho(Tr.Jacobian(), nor); + } + + double jcoef; + double wLM; + + const double w = ip.weight / Tr.Elem1->Weight(); + wLM = model->EvalDGWeight(w, *Tr.Elem1, eip); + jcoef = kappa * wLM * (nor*nor); + + for (int im = 0, i = 0; im < dim; ++im) + { + const double tj = jcoef * u_dir(im); + for (int idof = 0; idof < ndofs; ++idof, ++i) + { + elvect(i) += tj*shape(idof); + } + } + + } + +}; + +} // namespace mfem diff --git a/src/nlelast_solver.cpp b/src/nlelast_solver.cpp new file mode 100644 index 00000000..cdabc58b --- /dev/null +++ b/src/nlelast_solver.cpp @@ -0,0 +1,715 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include "nlelast_solver.hpp" +#include "input_parser.hpp" +#include "linalg_utils.hpp" +#include "etc.hpp" + +using namespace std; +using namespace mfem; + + +/* NLElastOperator */ + +NLElastOperator::NLElastOperator(const int height_, const int width_, Array &hs_, InterfaceForm *nl_itf_, + Array &u_offsets_, const bool direct_solve_) + : Operator(height_, width_), hs(hs_), nl_itf(nl_itf_), + u_offsets(u_offsets_), direct_solve(direct_solve_) +{ + + // TODO: this needs to be changed for parallel implementation. + sys_glob_size = height; + sys_row_starts[0] = 0; + sys_row_starts[1] = height; + + Hop = new BlockOperator(u_offsets); + for (int m = 0; m < hs_.Size(); m++) + { + assert(hs_[m]); + Hop->SetDiagonalBlock(m, hs_[m]); + } + + hs_mats.SetSize(hs.Size(), hs.Size()); + hs_mats = NULL; +} + +NLElastOperator::~NLElastOperator() +{ + delete Hop; + delete system_jac; + delete hs_jac; + delete uu_mono; + delete mono_jac; + delete jac_hypre; +} + +void NLElastOperator::Mult(const Vector &x, Vector &y) const +{ + y = 0.0; + //cout<<"in mult"<Mult(x, y); + //cout<<"out mult"<InterfaceAddMult(x_u, y_u); // TODO: Add this when interface integrator is there +} + +Operator& NLElastOperator::GetGradient(const Vector &x) const +{ + //cout<<"in grad"<(x), u_offsets[i], u_offsets[i+1] - u_offsets[i]); + //cout<<"before getgrad"<(&hs[i]->GetGradient(x_u)); + //cout<<"after getgrad"<(x), 0, u_offsets.Last()); + //if (nl_itf) nl_itf->InterfaceGetGradient(x_u, hs_mats); // TODO: Enable when we have interface integrator + + for (int i = 0; i < hs.Size(); i++) + for (int j = 0; j < hs.Size(); j++) + { + hs_mats(i, j)->Finalize(); + hs_jac->SetBlock(i, j, hs_mats(i, j)); + } + + mono_jac = hs_jac->CreateMonolithic(); + //cout<<"out grad"<("discretization/interface/kappa", (order + 1) * (order + 1)); + + var_names = GetVariableNames(); + num_var = var_names.size(); + + model = _model; + + // solution dimension is determined by initialization. + udim = dim; + vdim.SetSize(num_var); + vdim = dim; + + // Set up FE collection/spaces. + fec.SetSize(num_var); + if (full_dg) + { + fec = new DG_FECollection(order, dim, BasisType::GaussLobatto); + } + else + { + fec = new H1_FECollection(order, dim); + } + + fes.SetSize(numSub); + for (int m = 0; m < numSub; m++) + { + fes[m] = new FiniteElementSpace(meshes[m], fec[0], udim); + } +} + +NLElastSolver::~NLElastSolver() +{ + delete a_itf; + + DeletePointers(bs); + DeletePointers(as); + DeletePointers(lambda_c); + DeletePointers(mu_c); + + delete globalMat_mono; + delete globalMat; + delete globalMat_hypre; + delete mumps; + delete init_x; + +} + +void NLElastSolver::SetupIC(std::function F) +{ + init_x = new VectorFunctionCoefficient(dim, F); + for (int m = 0; m < numSub; m++) + { + assert(us[m]); + us[m]->ProjectCoefficient(*init_x); + + /* for (size_t i = 0; i < us[m]->Size(); i++) + { + cout<<"us[m]->Elem(i) is: "<Elem(i)<GetTrueVSize(); + num_vdofs[i] = fes[i]->GetTrueVSize(); + for (int d = 0; d < udim; d++) + { + block_offsets[d + i * udim + 1] = fes[i]->GetNDofs(); + } + } + block_offsets.PartialSum(); + var_offsets.PartialSum(); + domain_offsets = var_offsets; + + SetupBCVariables(); + + // Set up solution/rhs variables/ + U = new BlockVector(var_offsets); + RHS = new BlockVector(var_offsets); + /* + Note: for compatibility with ROM, it's better to split with domain_offsets. + For vector-component operations, can set up a view BlockVector like below: + + BlockVector *U_blocks = new BlockVector(U->GetData(), block_offsets); + + U_blocks does not own the data. + These are system-specific, therefore not defining it now. + */ + us.SetSize(numSub); + for (int m = 0; m < numSub; m++) + { + us[m] = new GridFunction(fes[m], U->GetBlock(m), 0); + + // BC's are weakly constrained and there is no essential dofs. + // Does this make any difference? + meshes[m]->GetNodes(*us[m]); + us[m]->SetTrueVector(); + //PrintVector(U->GetBlock(m), "ub.txt"); + } + if (use_rom) + MultiBlockSolver::InitROMHandler(); +} + +void NLElastSolver::BuildOperators() +{ + BuildRHSOperators(); + BuildDomainOperators(); +} + +bool NLElastSolver::BCExistsOnBdr(const int &global_battr_idx) +{ + assert((global_battr_idx >= 0) && (global_battr_idx < global_bdr_attributes.Size())); + assert(bdr_coeffs.Size() == global_bdr_attributes.Size()); + return (bdr_coeffs[global_battr_idx]); +} + +void NLElastSolver::BuildRHSOperators() +{ + bs.SetSize(numSub); + for (int m = 0; m < numSub; m++) + { + bs[m] = new LinearForm(fes[m], RHS->GetBlock(m).GetData()); + for (int r = 0; r < rhs_coeffs.Size(); r++) + { + bs[m]->AddDomainIntegrator(new VectorDomainLFIntegrator(*rhs_coeffs[r])); + } + } +} + +void NLElastSolver::SetupRHSBCOperators() +{ + assert(bs.Size() == numSub); + for (int m = 0; m < numSub; m++) + { + assert(bs[m]); + for (int b = 0; b < global_bdr_attributes.Size(); b++) + { + int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); + if (idx < 0) + continue; + + if (!BCExistsOnBdr(b)) + continue; + + switch (bdr_type[b]) + { + case BoundaryType::DIRICHLET: + //bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator( + //*bdr_coeffs[b], *lambda_c[m], *mu_c[m], alpha, kappa), *bdr_markers[b]); + bs[m]->AddBdrFaceIntegrator(new DGHyperelasticDirichletLFIntegrator( + *bdr_coeffs[b], model, alpha, kappa), *bdr_markers[b]); + break; + case BoundaryType::NEUMANN: + bs[m]->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(*bdr_coeffs[b]), *bdr_markers[b]); + break; + default: + printf("NLElastSolver::SetupRHSBCOperators - "); + printf("boundary attribute %d has a non-zero function, but does not have boundary type. will not be enforced.", + global_bdr_attributes[b]); + break; + } + } + } +} + +void NLElastSolver::BuildDomainOperators() +{ + // SanityCheckOnCoeffs(); + as.SetSize(numSub); + + for (int m = 0; m < numSub; m++) + { + as[m] = new NonlinearForm(fes[m]); + as[m]->AddDomainIntegrator(new HyperelasticNLFIntegratorHR(model)); + + if (full_dg) + { + as[m]->AddInteriorFaceIntegrator( + new DGHyperelasticNLFIntegrator(model, alpha, kappa)); + } + } + + a_itf = new InterfaceForm(meshes, fes, topol_handler); // TODO: Is this reasonable? + //a_itf->AddIntefaceIntegrator(new InterfaceDGElasticityIntegrator(lambda_c[0], mu_c[0], alpha, kappa)); +} + +void NLElastSolver::Assemble() +{ + AssembleRHS(); + //AssembleOperator(); +} + +void NLElastSolver::AssembleRHS() +{ + // SanityCheckOnCoeffs(); + MFEM_ASSERT(bs.Size() == numSub, "LinearForm bs != numSub.\n"); + for (int m = 0; m < numSub; m++) + { + MFEM_ASSERT(bs[m], "LinearForm or BilinearForm pointer of a subdomain is not associated!\n"); + bs[m]->Assemble(); + } + + for (int m = 0; m < numSub; m++) + // Do we really need SyncAliasMemory? + bs[m]->SyncAliasMemory(*RHS); // Synchronize with block vector RHS. What is different from SyncMemory? +} + +void NLElastSolver::AssembleOperator() +{ + "AssembleOperator is not needed for NLElastSolver!\n"; +} + +void NLElastSolver::AssembleInterfaceMatrices() +{ + assert(a_itf); + a_itf->AssembleInterfaceMatrices(mats); +} + +bool NLElastSolver::Solve() +{ + int maxIter = config.GetOption("solver/max_iter", 100); + double rtol = config.GetOption("solver/relative_tolerance", 1.e-10); + double atol = config.GetOption("solver/absolute_tolerance", 1.e-10); + int print_level = config.GetOption("solver/print_level", 0); + + int jac_maxIter = config.GetOption("solver/jacobian/max_iter", 10000); + double jac_rtol = config.GetOption("solver/jacobian/relative_tolerance", 1.e-10); + double jac_atol = config.GetOption("solver/jacobian/absolute_tolerance", 1.e-10); + int jac_print_level = config.GetOption("solver/jacobian/print_level", -1); + + bool lbfgs = config.GetOption("solver/use_lbfgs", false); + bool use_restart = config.GetOption("solver/use_restart", false); + std::string restart_file; + if (use_restart) + restart_file = config.GetRequiredOption("solver/restart_file"); + + const int hw = U->Size(); + NLElastOperator oper(hw, hw, as, a_itf, var_offsets, direct_solve); + + // Test mult + /* for (int i = 0; i < var_offsets.Size(); i++) + { + cout<<"var_offsets(i) is: "<GetBlock(0)); + Vector _U0(_U); + Vector _RHS(RHS->GetBlock(0)); + Vector _res(_U); + _res = 0.0; + cout<<"_U.Norml2() is: "<<_U.Norml2()< 0.0000001) + { + cout<<"ohnon "<< abs(_res(i) - _RHS(i))< 0.0000001) + { + cout<<"ohnon "<< abs(_res(i) - _RHS(i))<SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC); + mumps->SetPrintLevel(jac_print_level); + J_solver = mumps; + } + else + { + J_gmres = new GMRESSolver; + J_gmres->SetAbsTol(jac_atol); + J_gmres->SetRelTol(jac_rtol); + J_gmres->SetMaxIter(jac_maxIter); + J_gmres->SetPrintLevel(jac_print_level); + J_solver = J_gmres; + } + //cout<<"2 "<SetSolver(*J_solver); + newton_solver->SetOperator(oper); + newton_solver->SetPrintLevel(print_level); // print Newton iterations + newton_solver->SetRelTol(rtol); + newton_solver->SetAbsTol(atol); + newton_solver->SetMaxIter(maxIter); + cout<<"4 "<Mult(*RHS, *U); + //cout<<"5 "<GetConverged(); + //cout<<"6 "< F, const int battr) +{ + assert(bdr_coeffs.Size() > 0); + + if (battr > 0) + { + int idx = global_bdr_attributes.Find(battr); + if (idx < 0) + { + std::string msg = "battr " + std::to_string(battr) + " is not in global boundary attributes. skipping this boundary condition.\n"; + mfem_warning(msg.c_str()); + return; + } + bdr_coeffs[idx] = new VectorFunctionCoefficient(dim, F); + } + else + for (int k = 0; k < bdr_coeffs.Size(); k++) + bdr_coeffs[k] = new VectorFunctionCoefficient(dim, F); +} + +void NLElastSolver::AddRHSFunction(std::function F) +{ + rhs_coeffs.Append(new VectorFunctionCoefficient(dim, F)); +} + +void NLElastSolver::SetupBCOperators() +{ + SetupRHSBCOperators(); + SetupDomainBCOperators(); +} + +void NLElastSolver::SetupDomainBCOperators() +{ + MFEM_ASSERT(as.Size() == numSub, "NonlinearForm bs != numSub.\n"); + if (full_dg) + { + for (int m = 0; m < numSub; m++) + { + for (int b = 0; b < global_bdr_attributes.Size(); b++) + { + int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); + if (idx < 0) + continue; + if (!BCExistsOnBdr(b)) + continue; + + if (bdr_type[b] == BoundaryType::DIRICHLET) + //as[m]->AddBdrFaceIntegrator(new DGElasticityIntegrator( + // *(lambda_c[m]), *(mu_c[m]), alpha, kappa), *(bdr_markers[b])); + //cout<<"*lambda_c[m] is: "<lambda<mu<AddBdrFaceIntegrator( + new DGHyperelasticNLFIntegrator(model, alpha, kappa), *(bdr_markers[b])); + } + } + } +} + +void NLElastSolver::SetParameterizedProblem(ParameterizedProblem *problem) +{ + /* set up boundary types */ + MultiBlockSolver::SetParameterizedProblem(problem); + + // Set materials + lambda_c.SetSize(numSub); + lambda_c = NULL; + + mu_c.SetSize(numSub); + mu_c = NULL; + + Vector _x(1); + + for (size_t i = 0; i < numSub; i++) + { + double lambda_i = (problem->general_scalar_ptr[0])(_x); + lambda_c[i] = new ConstantCoefficient(lambda_i); + + double mu_i = (problem->general_scalar_ptr[1])(_x); + mu_c[i] = new ConstantCoefficient(mu_i); + } + + // Set BCs, the switch on BC type is done inside SetupRHSBCOperators + for (int b = 0; b < problem->battr.Size(); b++) + { + /* Dirichlet bc requires a function specified, even for zero. */ + if (problem->bdr_type[b] == BoundaryType::DIRICHLET) + assert(problem->vector_bdr_ptr[b]); + + /* Neumann bc does not require a function specified for zero */ + if (problem->vector_bdr_ptr[b]) + AddBCFunction(*(problem->vector_bdr_ptr[b]), problem->battr[b]); + } + + // Set RHS + if (problem->vector_rhs_ptr != NULL){ + AddRHSFunction(*(problem->vector_rhs_ptr)); + } + + // Add initial condition + if (problem->general_vector_ptr[0] != NULL) + { + SetupIC(*(problem->general_vector_ptr[0])); + } +} + +void NLElastSolver::ProjectOperatorOnReducedBasis() +{ + Array2D tmp(mats.NumRows(), mats.NumCols()); + for (int i = 0; i < tmp.NumRows(); i++) + for (int j = 0; j < tmp.NumCols(); j++) + tmp(i, j) = mats(i, j); + + rom_handler->ProjectOperatorOnReducedBasis(tmp); +} + +// Component-wise assembly +void NLElastSolver::BuildCompROMElement(Array &fes_comp) +{ + assert(train_mode == UNIVERSAL); + assert(rom_handler->BasisLoaded()); + + const int num_comp = fes_comp.Size(); + assert(comp_mats.Size() == num_comp); + + for (int c = 0; c < num_comp; c++) + { + Mesh *comp = topol_handler->GetComponentMesh(c); + BilinearForm a_comp(fes_comp[c]); + + a_comp.AddDomainIntegrator(new ElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]))); + if (full_dg) + a_comp.AddInteriorFaceIntegrator(new DGElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]), alpha, kappa)); + + a_comp.Assemble(); + a_comp.Finalize(); + + // Elasticity equation has only one solution variable. + comp_mats[c]->SetSize(1, 1); + (*comp_mats[c])(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); + } +} + +void NLElastSolver::BuildBdrROMElement(Array &fes_comp) +{ + assert(train_mode == UNIVERSAL); + assert(rom_handler->BasisLoaded()); + + const int num_comp = fes_comp.Size(); + assert(bdr_mats.Size() == num_comp); + + for (int c = 0; c < num_comp; c++) + { + Mesh *comp = topol_handler->GetComponentMesh(c); + assert(bdr_mats[c]->Size() == comp->bdr_attributes.Size()); + + MatrixBlocks *bdr_mat; + for (int b = 0; b < comp->bdr_attributes.Size(); b++) + { + Array bdr_marker(comp->bdr_attributes.Max()); + bdr_marker = 0; + bdr_marker[comp->bdr_attributes[b] - 1] = 1; + BilinearForm a_comp(fes_comp[c]); + a_comp.AddBdrFaceIntegrator(new DGElasticityIntegrator(*(lambda_c[c]), *(mu_c[c]), alpha, kappa), bdr_marker); + + a_comp.Assemble(); + a_comp.Finalize(); + + bdr_mat = (*bdr_mats[c])[b]; + bdr_mat->SetSize(1, 1); + (*bdr_mat)(0, 0) = rom_handler->ProjectToRefBasis(c, c, &(a_comp.SpMat())); + } + } +} + +void NLElastSolver::BuildInterfaceROMElement(Array &fes_comp) +{ + assert(topol_mode == TopologyHandlerMode::COMPONENT); + assert(train_mode == UNIVERSAL); + assert(rom_handler->BasisLoaded()); + + const int num_ref_ports = topol_handler->GetNumRefPorts(); + assert(port_mats.Size() == num_ref_ports); + for (int p = 0; p < num_ref_ports; p++) + { + assert(port_mats[p]->nrows == 2); + assert(port_mats[p]->ncols == 2); + + int c1, c2; + topol_handler->GetComponentPair(p, c1, c2); + + Array c_idx(2); + c_idx[0] = c1; + c_idx[1] = c2; + + Array2D spmats(2,2); + spmats = NULL; + + // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. + // Need to use two copied instances. + a_itf->AssembleInterfaceMatrixAtPort(p, fes_comp, spmats); + + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + (*port_mats[p])(i, j) = rom_handler->ProjectToRefBasis(c_idx[i], c_idx[j], spmats(i,j)); + + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) delete spmats(i, j); + } // for (int p = 0; p < num_ref_ports; p++) +} +//void NLElastSolver::SanityCheckOnCoeffs() { "NLElastSolver::SanityCheckOnCoeffs is not implemented yet!\n"; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 2f2ec3fd..44045a8e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -14,8 +14,10 @@ add_executable(poisson_dd_mms poisson_dd_mms.cpp $ add_executable(stokes_dd_mms stokes_dd_mms.cpp $ $) add_executable(steady_ns_dd_mms steady_ns_dd_mms.cpp $ $) add_executable(linelast_dd_mms linelast_dd_mms.cpp $ $) +add_executable(nlelast_dd_mms nlelast_dd_mms.cpp $ $) add_executable(advdiff_dd_mms advdiff_dd_mms.cpp $ $) add_executable(dg_integ_mms dg_integ_mms.cpp $ $) +add_executable(neohookean_test neohookean_test.cpp $ $) file(COPY inputs/dd_mms.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) file(COPY meshes/dd_mms.mesh DESTINATION ${CMAKE_BINARY_DIR}/test/meshes) file(COPY inputs/dd_mms.component.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs) @@ -71,6 +73,9 @@ add_executable(test_rom_nonlinearform test_rom_nonlinearform.cpp $) +add_executable(nlelast_test nlelast_test.cpp $) + + function(add_test_dir TEST_DIR) # helper function that adds all executables in a given directory as tests for ctest get_property(target_names DIRECTORY ${TEST_DIR} PROPERTY BUILDSYSTEM_TARGETS) diff --git a/test/linelast_dd_mms.cpp b/test/linelast_dd_mms.cpp index d27f9152..d8a537aa 100644 --- a/test/linelast_dd_mms.cpp +++ b/test/linelast_dd_mms.cpp @@ -7,7 +7,7 @@ using namespace std; using namespace mfem; -using namespace mms::linelast; +using namespace mms::nlelast; /** * Simple smoke test to make sure Google Test is properly linked @@ -17,17 +17,34 @@ TEST(GoogleTestFramework, GoogleTestFrameworkFound) SUCCEED(); } +/* TEST(DDSerialTest, Test_convergence_DG_integratorwise) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["mesh"]["filename"] = "meshes/test.1x1.mesh"; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "none"; + bool nonlinear = true; + CheckConvergenceIntegratorwise(); + + return; +} */ + TEST(DDSerialTest, Test_convergence_DG) { config = InputParser("inputs/dd_mms.yml"); + //config.dict_["mesh"]["filename"] = "meshes/test.1x1.mesh"; config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; + config.dict_["solver"]["direct_solve"] = true; config.dict_["discretization"]["full-discrete-galerkin"] = true; config.dict_["domain-decomposition"]["type"] = "none"; - CheckConvergence(); + bool nonlinear = false; + CheckConvergence(nonlinear); return; } +/* TEST(DDSerialTest, Test_direct_solver_DG) { config = InputParser("inputs/dd_mms.yml"); @@ -35,12 +52,29 @@ TEST(DDSerialTest, Test_direct_solver_DG) config.dict_["solver"]["direct_solve"] = true; config.dict_["discretization"]["full-discrete-galerkin"] = true; config.dict_["domain-decomposition"]["type"] = "none"; - CheckConvergence(); + bool nonlinear = true; + CheckConvergence(nonlinear); return; -} +} */ -TEST(DDSerialTest, Test_convergence_DG_DD) +/* TEST(DDSerialTest, CompareSolvers) +{ + config = InputParser("inputs/dd_mms.yml"); + //config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; + config.dict_["mesh"]["filename"] = "meshes/test.1x1.mesh"; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "none"; + bool nonlinear = false; + config.dict_["discretization"]["interface/alpha"] = 0.0; + + CompareLinMat(); + return; +} */ + + +/* TEST(DDSerialTest, Test_convergence_DG_DD) { config = InputParser("inputs/dd_mms.yml"); config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; @@ -49,9 +83,9 @@ TEST(DDSerialTest, Test_convergence_DG_DD) CheckConvergence(); return; -} +} */ -TEST(DDSerialTest, Test_direct_solver_DG_DD) +/* TEST(DDSerialTest, Test_direct_solver_DG_DD) { config = InputParser("inputs/dd_mms.yml"); config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; @@ -61,7 +95,7 @@ TEST(DDSerialTest, Test_direct_solver_DG_DD) CheckConvergence(); return; -} +} */ int main(int argc, char *argv[]) { diff --git a/test/mms_suite.cpp b/test/mms_suite.cpp index 5d95222f..61f45535 100644 --- a/test/mms_suite.cpp +++ b/test/mms_suite.cpp @@ -597,6 +597,480 @@ namespace linelast } // namespace linelast +namespace nlelast +{ + void ExactSolutionLinear(const Vector &x, Vector &u) + { + u = 0.0; + for (size_t i = 0; i < dim; i++) + { + u(i) = pow(x(i), 3.0); + } + } + + void ExactRHSLinear(const Vector &x, Vector &u) + { + u = 0.0; + for (size_t i = 0; i < dim; i++) + { + u(i) = 6.0 * x(i) * (lambda + 2.0 * mu); + } + u *= -1.0; + } + + void ExactSolutionNeoHooke(const Vector &x, Vector &u) + { + u = 0.0; + //assert(dim == 2); + assert(x.Size() == 2); + u(0) = pow(x(0), 2.0) + x(0); + u(1) = pow(x(1), 2.0) + x(1); + } + + void ExactSolutionNeoHookeBC(const Vector &x, Vector &u) + { + u = 0.0; + //assert(dim == 2); + assert(x.Size() == 2); + u(0) = pow(x(0), 2.0); + u(1) = pow(x(1), 2.0); + } + + void SimpleExactSolutionNeoHooke(const Vector &X, Vector &U) + { + int dim = 2; + int dof = U.Size()/dim; + U = 0.0; + for (size_t i = 0; i < U.Size()/dim; i++) + { + U(i) = pow(X(i), 2.0) + X(i); + U(dof + i) = pow(X(dof + i), 2.0) + X(dof + i); + } + } + + /* void ExactRHSNeoHooke(const Vector &x, Vector &u) + { + u = 0.0; + assert(dim == 2); + + const double x_1 = x(0); + const double x_2 = x(1); + + u(0) = (128.0*K + 128.0*mu + 1024.0*K*x_1 + 1024.0*K*x_2 + 128.0*mu*x_1 + 384.0*mu*x_2 + 3072.0*K*(pow(x_1,2)) + 8192.0*K*x_1*x_2 + 3072.0*K*(pow(x_2,2)) + 128.0*mu*(pow(x_1,2)) + 384.0*mu*(pow(x_2,2)) + 4096.0*K*(pow(x_1,3)) + 24576.0*K*(pow(x_1,2))*x_2 + 24576.0*K*x_1*(pow(x_2,2)) + 4096.0*K*(pow(x_2,3)) + 2048.0*K*(pow(x_1,4)) + 32768.0*K*(pow(x_1,3))*x_2 + 73728.0*K*(pow(x_1,2))*(pow(x_2,2)) + 32768.0*K*x_1*(pow(x_2,3)) + 2048.0*K*(pow(x_2,4)) + 16384.0*K*(pow(x_1,4))*x_2 + 98304.0*K*(pow(x_1,3))*(pow(x_2,2)) + 98304.0*K*(pow(x_1,2))*(pow(x_2,3)) + 16384.0*K*x_1*(pow(x_2,4)) + 49152.0*K*(pow(x_1,4))*(pow(x_2,2)) + 131072.0*K*(pow(x_1,3))*(pow(x_2,3)) + 49152.0*K*(pow(x_1,2))*(pow(x_2,4)) + 65536.0*K*(pow(x_1,4))*(pow(x_2,3)) + 65536.0*K*(pow(x_1,3))*(pow(x_2,4)) + 32768.0*K*(pow(x_1,4))*(pow(x_2,4))) / (64.0 * (pow((1 + 2*x_1), 4))*(pow((1 + 2*x_2),2))); + + u(1) = (128.0*K + 128.0*mu + 1024.0*K*x_1 + 1024.0*K*x_2 + 384.0*mu*x_1 + 128.0*mu*x_2 + 3072.0*K*(pow(x_1,2)) + 8192.0*K*x_1*x_2 + 3072.0*K*(pow(x_2,2)) + 384.0*mu*(pow(x_1,2)) + 128.0*mu*(pow(x_2,2)) + 4096.0*K*(pow(x_1,3)) + 24576.0*K*(pow(x_1,2))*x_2 + 24576.0*K*x_1*(pow(x_2,2)) + 4096.0*K*(pow(x_2,3)) + 2048.0*K*(pow(x_1,4)) + 32768.0*K*(pow(x_1,3))*x_2 + 73728.0*K*(pow(x_1,2))*(pow(x_2,2)) + 32768.0*K*x_1*(pow(x_2,3)) + 2048.0*K*(pow(x_2,4)) + 16384.0*K*(pow(x_1,4))*x_2 + 98304.0*K*(pow(x_1,3))*(pow(x_2,2)) + 98304.0*K*(pow(x_1,2))*(pow(x_2,3)) + 16384.0*K*x_1*(pow(x_2,4)) + 49152.0*K*(pow(x_1,4))*(pow(x_2,2)) + 131072.0*K*(pow(x_1,3))*(pow(x_2,3)) + 49152.0*K*(pow(x_1,2))*(pow(x_2,4)) + 65536.0*K*(pow(x_1,4))*(pow(x_2,3)) + 65536.0*K*(pow(x_1,3))*(pow(x_2,4)) + 32768.0*K*(pow(x_1,4))*(pow(x_2,4))) / (64.0 * (pow((1 + 2*x_1),2))*(pow((1 + 2*x_2),4))); + //u *= -1.0; + } */ + + void SimpleExactRHSNeoHooke(const Vector &x, Vector &u) + { + u = 0.0; + assert(dim == 2); + assert(mu == 0.0); + u(0) = 2 * K * pow(1.0 + 2.0 * x(1), 2.0); + u(1) = 2 * K * pow(1.0 + 2.0 * x(0), 2.0); + u *= -1.0; + } + + void NullSolution(const Vector &x, Vector &u) + { + u = 0.0; + //u -= x; + } + + void NullDefSolution(const Vector &x, Vector &u) + { + u = 0.0; + u = x; + //u = 1.0; + } + + void cantileverf(const Vector &x, Vector &u) + { + u = 0.0; + u(0) = -0.10; + } + + void cantileverfu(const Vector &x, Vector &u) + { + u = 0.0; + u(0) = 0.10; + } + + NLElastSolver *SolveWithRefinement(const int num_refinement, const bool nonlinear) + { + config.dict_["mesh"]["uniform_refinement"] = num_refinement; + DGHyperelasticModel *model = NULL; + + if (nonlinear) + { + model = new NeoHookeanHypModel(mu, K); + + } + else + { + model = new LinElastMaterialModel(mu, lambda); + } + + NLElastSolver *test = new NLElastSolver(model); + + + dim = test->GetDim(); + + test->InitVariables(); + test->InitVisualization(); + if (nonlinear) + { + test->AddBCFunction(ExactSolutionNeoHooke); + //test->AddBCFunction(ExactSolutionNeoHookeBC, 1); + //test->AddBCFunction(ExactSolutionNeoHookeBC, 2); + //test->AddBCFunction(ExactSolutionNeoHookeBC, 3); + //test->AddBCFunction(cantileverf, 1); + //test->AddBCFunction(cantileverfu, 2); + //test->AddBCFunction(NullSolution, 3); + test->AddRHSFunction(SimpleExactRHSNeoHooke); + //test->SetupIC(NullDefSolution); + test->SetBdrType(BoundaryType::DIRICHLET); + + } + else + { + test->AddBCFunction(ExactSolutionLinear); + /* test->AddBCFunction(ExactSolutionLinear, 2); + test->AddBCFunction(ExactSolutionLinear, 3); */ + test->AddRHSFunction(ExactRHSLinear); + test->SetupIC(ExactSolutionLinear); + } + //test->AddBCFunction(NullSolution, 1); + //test->AddBCFunction(NullSolution, 2); + //test->AddBCFunction(NullSolution, 3); + //test->AddBCFunction(cantileverfu, 2); + //test->AddBCFunction(NullSolution, 3); + test->SetBdrType(BoundaryType::DIRICHLET); + //test->SetBdrType(BoundaryType::NEUMANN,1); + //test->SetBdrType(BoundaryType::NEUMANN,2); + //test->SetupIC(ExactSolutionNeoHooke); + + test->BuildOperators(); + + test->SetupBCOperators(); + + test->Assemble(); + + test->Solve(); + + return test; + } + + + void CompareLinMat() + { + int num_refine = config.GetOption("manufactured_solution/number_of_refinement", 3); + int base_refine = config.GetOption("manufactured_solution/baseline_refinement", 0); + + // Compare with exact solution + config.dict_["mesh"]["uniform_refinement"] = 0; + LinElastMaterialModel* model = new LinElastMaterialModel(mu, lambda); + NLElastSolver *test1 = new NLElastSolver(model); + + LinElastSolver *test2 = new LinElastSolver(mu, lambda); + dim = test2->GetDim(); + assert(dim == 2); + test2->InitVariables(); + test2->InitVisualization(); + test2->AddBCFunction(ExactSolutionLinear, 1); + test2->AddBCFunction(ExactSolutionLinear, 2); + test2->AddBCFunction(ExactSolutionLinear, 3); + test2->SetBdrType(BoundaryType::DIRICHLET); + test2->AddRHSFunction(ExactRHSLinear); + test2->BuildOperators(); + test2->SetupBCOperators(); + test2->Assemble(); + test2->Solve(); + + dim = test1->GetDim(); + assert(dim == 2); + test1->InitVariables(); + test1->InitVisualization(); + test1->AddBCFunction(ExactSolutionLinear, 1); + test1->AddBCFunction(ExactSolutionLinear, 2); + test1->AddBCFunction(ExactSolutionLinear, 3); + test1->SetBdrType(BoundaryType::DIRICHLET); + test1->AddRHSFunction(ExactRHSLinear); + test1->BuildOperators(); + test1->SetupBCOperators(); + test1->Assemble(); + test1->Solve(); + + return; + } + + void CheckConvergence(const bool nonlinear) + { + int num_refine = config.GetOption("manufactured_solution/number_of_refinement", 3); + int base_refine = config.GetOption("manufactured_solution/baseline_refinement", 0); + + // Compare with exact solution + int dim = 2; // only check two dimensions + //const double mu = 0.0; + //const double K = 1.0; + VectorFunctionCoefficient* exact_sol; + if (nonlinear) + { + exact_sol = new VectorFunctionCoefficient(dim, ExactSolutionNeoHooke); + } + else + { + exact_sol = new VectorFunctionCoefficient(dim, ExactSolutionLinear); + } + + printf("Num. Elem.\tRelative Error\tConvergence Rate\tNorm\n"); + + Vector conv_rate(num_refine); + conv_rate = 0.0; + double error1 = 0.0; + base_refine = 0; + num_refine = 4; + for (int r = base_refine; r < num_refine; r++) + { + NLElastSolver *test = SolveWithRefinement(r, nonlinear); + + int order = test->GetDiscretizationOrder(); + cout<<"order is: "<GetNumSubdomains(); k++) + { + Mesh *mk = test->GetMesh(k); + norm += pow(ComputeLpNorm(2.0, *exact_sol, *mk, irs), 2); + numEl += mk->GetNE(); + } + norm = sqrt(norm); + double error = 0.0; + for (int k = 0; k < test->GetNumSubdomains(); k++) + { + GridFunction *uk = test->GetGridFunction(k); + error += pow(uk->ComputeLpError(2, *exact_sol), 2); + } + error = sqrt(error); + error /= norm; + + if (r > base_refine) + { + conv_rate(r) = error1 / error; + } + + printf("%d\t%.15E\t%.15E\t%.15E\n", numEl, error, conv_rate(r), norm); + + // reported convergence rate + if (r > base_refine) + //EXPECT_TRUE(conv_rate(r) > pow(2.0, order + 1) - 0.5); + EXPECT_TRUE(conv_rate(r) > pow(2.0, order + 1) - 0.8); + + error1 = error; + } + + return; + } + + void test_fn(const Vector &x, Vector &u) +{ + double xi(x(0)); + double yi(x(1)); + + assert(x.Size() == 2); + + u(0) = sin(pi * x(0)); + u(1) = sin(pi * x(1)); + return; +} + + double EvalWithRefinement(const int num_refinement, int &order_out) +{ + // 1. Parse command-line options. + std::string mesh_file = config.GetRequiredOption("mesh/filename"); + bool use_dg = config.GetOption("discretization/full-discrete-galerkin", false); + int order = config.GetOption("discretization/order", 1); + order_out = order; + + Mesh *mesh = new Mesh(mesh_file.c_str(), 1, 1); + int dim = mesh->Dimension(); + //cout<<"dim is: "<UniformRefinement(); + } + + FiniteElementCollection *dg_coll(new DG_FECollection(order, dim)); + FiniteElementCollection *h1_coll(new H1_FECollection(order, dim)); + + FiniteElementSpace *fes; + //FiniteElementSpace fespace(mesh, &fe_coll, dim); + if (use_dg) + { + fes = new FiniteElementSpace(mesh, dg_coll, dim); + } + else + { + fes = new FiniteElementSpace(mesh, h1_coll, dim); + } + + // 12. Create the grid functions u and p. Compute the L2 error norms. + //cout<<"pre error1"<GetTrueVSize(); + if (test_integ == "domain") + { + assert(use_dg == true); + + Array p_ess_attr(mesh->bdr_attributes.Max()); + // this array of integer essentially acts as the array of boolean: + // If value is 0, then it is not Dirichlet. + // If value is 1, then it is Dirichlet. + p_ess_attr = 1; + double kappa = -1.0; + //kappa = 0.0; + LinearForm *gform = new LinearForm(fes); + VectorFunctionCoefficient ud(dim, ExactSolutionNeoHooke); + + gform->AddBdrFaceIntegrator(new DGHyperelasticDirichletLFIntegrator( + ud, &model2, 0.0, kappa), p_ess_attr); + gform->Assemble(); + + NonlinearForm *nlform = new NonlinearForm(fes); + nlform->AddDomainIntegrator(new HyperelasticNLFIntegratorHR(&model2)); + nlform->AddBdrFaceIntegrator( new DGHyperelasticNLFIntegrator(&model2, 0.0, kappa),p_ess_attr); + + GridFunction x_ref(fes); + mesh->GetNodes(x_ref); + x.SetSize(ndofs); + x = x_ref.GetTrueVector(); + + y0.SetSize(ndofs); + y0 = 0.0; + SimpleExactSolutionNeoHooke(x, y0); + + y1.SetSize(ndofs); + y1 = 0.0; + nlform->Mult(y0, y1); //MFEM Neohookean + + for (size_t i = 0; i < y1.Size(); i++) + { + y1(i) -= gform->Elem(i); + } + + product = p * y1; + delete nlform; + delete gform; + } + else if (test_integ == "bc") + { + assert(use_dg == true); + Array p_ess_attr(mesh->bdr_attributes.Max()); + // this array of integer essentially acts as the array of boolean: + // If value is 0, then it is not Dirichlet. + // If value is 1, then it is Dirichlet. + p_ess_attr = 1; + //p_ess_attr[1] = 1; + LinearForm *gform = new LinearForm(fes); + VectorFunctionCoefficient ud(dim, ExactSolutionNeoHooke); + + gform->AddBdrFaceIntegrator(new DGHyperelasticDirichletLFIntegrator( + ud, &model2, 0.0, -1.0), p_ess_attr); + gform->Assemble(); + + NonlinearForm *nlform = new NonlinearForm(fes); + nlform->AddBdrFaceIntegrator(new DGHyperelasticNLFIntegrator(&model2, 0.0, -1.0), p_ess_attr); + y0.SetSize(ndofs); + y0 = 0.0; + SimpleExactSolutionNeoHooke(x, y0); + + y1.SetSize(ndofs); + y1 = 0.0; + nlform->Mult(y0, y1); //MFEM Neohookean + + //y1 -= gform; + product = y1.Norml2(); + delete gform; + + } + + // 17. Free the used memory. + delete fes; + delete dg_coll; + delete h1_coll; + delete mesh; + + return product; +} + +void CheckConvergenceIntegratorwise() +{ + int num_refine = config.GetOption("manufactured_solution/number_of_refinement", 5); + num_refine = 8; + //double product_ex = 26.0 * K / 3.0 * 1.5384588; // TODO: replace + string test_integ = "domain"; + double product_ex =0.0; + if (test_integ == "bc") + { + double wlm = K; + //product_ex = 4.0 * kappa * wlm * (pow(pi,2.0) - 4.0)/pow(pi,3.0); + //product_ex = 1.0; + product_ex = 4.0 * -1.0 * (pow(pi,2.0) - 4.0)/pow(pi,3.0); + + } + else if (test_integ == "domain") + { + product_ex = -(104.0 * K)/(3.0 * pi); // TODO: replace + } + + + printf("(p, n dot u_d)_ex = %.5E\n", product_ex); + + printf("Num. Refine.\tRel. Error\tConv Rate\tProduct\tProduct_ex\n"); + + Vector conv_rate(num_refine); + conv_rate = 0.0; + double error1 = 0.0; + for (int r = 0; r < num_refine; r++) + { + int order = -1; + double product = EvalWithRefinement(r, order); + + double error = abs(product - product_ex) / abs(product_ex); + + if (r > 0) + conv_rate(r) = error1 / error; + printf("%d\t%.5E\t%.5E\t%.5E\t%.5E\n", r, error, conv_rate(r), product, product_ex); + + // reported convergence rate + if (r > 0) + EXPECT_TRUE(conv_rate(r) > pow(2.0, order+1) - 0.1); + + error1 = error; + } + + return; +} + +} // namespace nlelast + namespace advdiff { diff --git a/test/mms_suite.hpp b/test/mms_suite.hpp index 0e49fc23..396d9fcc 100644 --- a/test/mms_suite.hpp +++ b/test/mms_suite.hpp @@ -10,6 +10,7 @@ #include "stokes_solver.hpp" #include "steady_ns_solver.hpp" #include "linelast_solver.hpp" +#include "nlelast_solver.hpp" #include "advdiff_solver.hpp" #include #include @@ -91,6 +92,25 @@ void CheckConvergence(); } // namespace linelast +namespace nlelast +{ + +static const double pi = 4.0 * atan(1.0); +//static const double mu = 3.14; +static const double mu = 1.0; +//static double K = 2.33; +static double K = 1.0; +static double lambda = K; +static int dim; +static void ExactSolution(const Vector & x, Vector & u); +static void ExactRHS(const Vector & x, Vector & u); +NLElastSolver *SolveWithRefinement(const int num_refinement); +void CheckConvergence(bool nonlinear); +void CheckConvergenceIntegratorwise(); +void CompareLinMat(); + +} // namespace nlelast + namespace advdiff { diff --git a/test/neohookean_test.cpp b/test/neohookean_test.cpp new file mode 100644 index 00000000..b308f027 --- /dev/null +++ b/test/neohookean_test.cpp @@ -0,0 +1,206 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include +#include +#include "mms_suite.hpp" +#include "nlelast_integ.hpp" + +using namespace std; +using namespace mfem; +namespace mfem +{ + void MFEMAssembleH(const DenseMatrix &J, const DenseMatrix &DS, + const double weight, DenseMatrix &A) + { + + /// set params + double g = 1.34; + double mu = 3.14; + double K = 4.13; + + int dof = DS.Height(), dim = DS.Width(); + + DenseMatrix Z(dim); + DenseMatrix G(dof, dim); + DenseMatrix C(dof, dim); + + double dJ = J.Det(); + double sJ = dJ/g; + double a = mu*pow(dJ, -2.0/dim); + double bc = a*(J*J)/dim; + double b = bc - K*sJ*(sJ - 1.0); + double c = 2.0*bc/dim + K*sJ*(2.0*sJ - 1.0); + + CalcAdjugateTranspose(J, Z); + Z *= (1.0/dJ); // Z = J^{-t} + + MultABt(DS, J, C); // C = DS J^t + MultABt(DS, Z, G); // G = DS J^{-1} + + a *= weight; + b *= weight; + c *= weight; + + // 1. + for (int i = 0; i < dof; i++) + for (int k = 0; k <= i; k++) + { + double s = 0.0; + for (int d = 0; d < dim; d++) + { + s += DS(i,d)*DS(k,d); + } + s *= a; + + for (int d = 0; d < dim; d++) + { + A(i+d*dof,k+d*dof) += s; + } + + if (k != i) + for (int d = 0; d < dim; d++) + { + A(k+d*dof,i+d*dof) += s; + } + } + + a *= (-2.0/dim); + + // 2. + for (int i = 0; i < dof; i++) + for (int j = 0; j < dim; j++) + for (int k = 0; k < dof; k++) + for (int l = 0; l < dim; l++) + { + A(i+j*dof,k+l*dof) += + a * (C(i,j)*G(k,l) + G(i,j)*C(k,l)) + + b*G(i,l)*G(k,j) + c*G(i,j)*G(k,l); + } + } + +void SimpleAssembleH(const DenseMatrix &J, const DenseMatrix &DS, + const double weight, DenseMatrix &A) +{ + /// set params + double g = 1.34; + double mu = 3.14; + double K = 4.13; + + int dof = DS.Height(), dim = DS.Width(); + + DenseMatrix Z(dim); + DenseMatrix G(dof, dim); + DenseMatrix C(dof, dim); + + double dJ = J.Det(); + double sJ = dJ/g; + double a = mu*pow(dJ, -2.0/dim); + double bc = a*(J*J)/dim; + double b = bc - K*sJ*(sJ - 1.0); + double c = 2.0*bc/dim + K*sJ*(2.0*sJ - 1.0); + + CalcAdjugateTranspose(J, Z); + Z *= (1.0/dJ); // Z = J^{-t} + + MultABt(DS, J, C); // C = DS J^t + MultABt(DS, Z, G); // G = DS J^{-1} + + a *= weight; + b *= weight; + c *= weight; + const double a2 = a * (-2.0/dim); + + for (size_t i = 0; i < dof; i++) + { + for (size_t j = 0; j < dim; j++) // Looping over each entry in residual + { + const int ij = j * dof + i; + + for (size_t m = 0; m < dof; m++) + for (size_t n = 0; n < dim; n++) // Looping over derivatives with respect to U + { + const int mn = n * dof + m; + double temp = 0.0; + for (size_t k = 0; k < dim; k++) + { + const int S_jk = k * dim + j; + //temp += Dmat(S_jk, mn) * w * gshape(i,k); + const double s1 = (j==n) ? a * DS(m,k) : 0.0; + const double s2 = a2 * (J(j,k)*G(m,n) + Z(j,k)*C(m,n)) + + b*Z(n,k)*G(m,j) + c*Z(j,k)*G(m,n); + + temp += DS(i,k)*(s1 + s2); + } + A(ij, mn) += temp; + } + } + } +} +} //namespace mfem + +/** + * Simple smoke test to make sure Google Test is properly linked + */ +TEST(GoogleTestFramework, GoogleTestFrameworkFound) { + SUCCEED(); +} + +// Check that SimpleAssembleH works +TEST(Assemble_H, Test_NLElast) +{ + int dim = 2; + int ndofs = 4; + + DenseMatrix J(dim); + DenseMatrix DS(ndofs, dim); + const double w = 1.2; + DenseMatrix A1(ndofs*dim, ndofs*dim); + DenseMatrix A2(ndofs*dim, ndofs*dim); + A1 = 0.0; + A2 = 0.0; + + double lower_bound = -1; + double upper_bound = 1; + + uniform_real_distribution unif(lower_bound, + upper_bound); + default_random_engine re; + + for (size_t i = 0; i < dim; i++) + for (size_t j = 0; j < dim; j++) + { + J(i,j) = unif(re); + } + + for (size_t i = 0; i < ndofs; i++) + for (size_t j = 0; j < dim; j++) + { + DS(i,j) = unif(re); + } + + MFEMAssembleH(J, DS, w, A1); + SimpleAssembleH(J, DS, w, A2); + + PrintMatrix(A1, "A1.txt"); + PrintMatrix(A2, "A2.txt"); + + cout << "MFEM Stiffness matrix norm: " << A1.FNorm() << endl; + cout << "ScaleupROM Stiffness matrix norm: " << A2.FNorm() << endl; + A1.Add(-1.0, A2); + double norm_diff = A1.FNorm(); + cout << "Stiffness matrix difference norm: " << norm_diff << endl; + + return; +} + + +int main(int argc, char *argv[]) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} \ No newline at end of file diff --git a/test/nlelast_dd_mms.cpp b/test/nlelast_dd_mms.cpp new file mode 100644 index 00000000..d8a537aa --- /dev/null +++ b/test/nlelast_dd_mms.cpp @@ -0,0 +1,107 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include +#include "mms_suite.hpp" + +using namespace std; +using namespace mfem; +using namespace mms::nlelast; + +/** + * Simple smoke test to make sure Google Test is properly linked + */ +TEST(GoogleTestFramework, GoogleTestFrameworkFound) +{ + SUCCEED(); +} + +/* TEST(DDSerialTest, Test_convergence_DG_integratorwise) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["mesh"]["filename"] = "meshes/test.1x1.mesh"; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "none"; + bool nonlinear = true; + CheckConvergenceIntegratorwise(); + + return; +} */ + +TEST(DDSerialTest, Test_convergence_DG) +{ + config = InputParser("inputs/dd_mms.yml"); + //config.dict_["mesh"]["filename"] = "meshes/test.1x1.mesh"; + config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "none"; + bool nonlinear = false; + CheckConvergence(nonlinear); + + return; +} + +/* +TEST(DDSerialTest, Test_direct_solver_DG) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "none"; + bool nonlinear = true; + CheckConvergence(nonlinear); + + return; +} */ + +/* TEST(DDSerialTest, CompareSolvers) +{ + config = InputParser("inputs/dd_mms.yml"); + //config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; + config.dict_["mesh"]["filename"] = "meshes/test.1x1.mesh"; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "none"; + bool nonlinear = false; + config.dict_["discretization"]["interface/alpha"] = 0.0; + + CompareLinMat(); + return; +} */ + + +/* TEST(DDSerialTest, Test_convergence_DG_DD) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "interior_penalty"; + CheckConvergence(); + + return; +} */ + +/* TEST(DDSerialTest, Test_direct_solver_DG_DD) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["mesh"]["filename"] = "../examples/linelast/meshes/beam-tri.mesh"; + config.dict_["solver"]["direct_solve"] = true; + config.dict_["discretization"]["full-discrete-galerkin"] = true; + config.dict_["domain-decomposition"]["type"] = "interior_penalty"; + CheckConvergence(); + + return; +} */ + +int main(int argc, char *argv[]) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} \ No newline at end of file diff --git a/test/nlelast_test.cpp b/test/nlelast_test.cpp new file mode 100644 index 00000000..615c9afe --- /dev/null +++ b/test/nlelast_test.cpp @@ -0,0 +1,655 @@ +// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: MIT + +#include +#include +#include "mms_suite.hpp" +#include "nlelast_integ.hpp" + +using namespace std; +using namespace mfem; + +namespace mfem +{ + class Test_DGElasticityIntegrator : public BilinearFormIntegrator + { + public: + Test_DGElasticityIntegrator(double alpha_, double kappa_) + : lambda(NULL), mu(NULL), alpha(alpha_), kappa(kappa_) { } + + Test_DGElasticityIntegrator(Coefficient &lambda_, Coefficient &mu_, + double alpha_, double kappa_) + : lambda(&lambda_), mu(&mu_), alpha(alpha_), kappa(kappa_) { } + + using BilinearFormIntegrator::AssembleFaceMatrix; + virtual void AssembleFaceMatrix(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Trans, + DenseMatrix &elmat); + + protected: + Coefficient *lambda, *mu; + double alpha, kappa; + + #ifndef MFEM_THREAD_SAFE + // values of all scalar basis functions for one component of u (which is a + // vector) at the integration point in the reference space + Vector shape1, shape2; + // values of derivatives of all scalar basis functions for one component + // of u (which is a vector) at the integration point in the reference space + DenseMatrix dshape1, dshape2; + // Adjugate of the Jacobian of the transformation: adjJ = det(J) J^{-1} + DenseMatrix adjJ; + // gradient of shape functions in the real (physical, not reference) + // coordinates, scaled by det(J): + // dshape_ps(jdof,jm) = sum_{t} adjJ(t,jm)*dshape(jdof,t) + DenseMatrix dshape1_ps, dshape2_ps; + Vector nor; // nor = |weight(J_face)| n + Vector nL1, nL2; // nL1 = (lambda1 * ip.weight / detJ1) nor + Vector nM1, nM2; // nM1 = (mu1 * ip.weight / detJ1) nor + Vector dshape1_dnM, dshape2_dnM; // dshape1_dnM = dshape1_ps . nM1 + // 'jmat' corresponds to the term: kappa + DenseMatrix jmat; + #endif + + static void AssembleBlock( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, + const double jmatcoef, const Vector &col_nL, const Vector &col_nM, + const Vector &row_shape, const Vector &col_shape, + const Vector &col_dshape_dnM, const DenseMatrix &col_dshape, + DenseMatrix &elmat, DenseMatrix &jmat); + }; + + void _AssembleBlock( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, + const double jmatcoef, const Vector &col_nL, const Vector &col_nM, + const Vector &row_shape, const Vector &col_shape, + const Vector &col_dshape_dnM, const DenseMatrix &col_dshape, + DenseMatrix &elmat, DenseMatrix &jmat) + { + for (int jm = 0, j = col_offset; jm < dim; ++jm) + { + for (int jdof = 0; jdof < col_ndofs; ++jdof, ++j) + { + const double t2 = col_dshape_dnM(jdof); + for (int im = 0, i = row_offset; im < dim; ++im) + { + const double t1 = col_dshape(jdof, jm) * col_nL(im); + const double t3 = col_dshape(jdof, im) * col_nM(jm); + const double tt = t1 + ((im == jm) ? t2 : 0.0) + t3; + for (int idof = 0; idof < row_ndofs; ++idof, ++i) + { + elmat(i, j) += row_shape(idof) * tt; + } + } + } + } + + if (jmatcoef == 0.0) { return; } + + for (int d = 0; d < dim; ++d) + { + const int jo = col_offset + d*col_ndofs; + const int io = row_offset + d*row_ndofs; + for (int jdof = 0, j = jo; jdof < col_ndofs; ++jdof, ++j) + { + const double sj = jmatcoef * col_shape(jdof); + for (int i = max(io,j), idof = i - io; idof < row_ndofs; ++idof, ++i) + { + jmat(i, j) += row_shape(idof) * sj; + } + } + } + }; + + void Test_DGElasticityIntegrator::AssembleFaceMatrix( + const FiniteElement &el1, const FiniteElement &el2, + FaceElementTransformations &Trans, DenseMatrix &elmat) + { + #ifdef MFEM_THREAD_SAFE + // For descriptions of these variables, see the class declaration. + Vector shape1, shape2; + DenseMatrix dshape1, dshape2; + DenseMatrix adjJ; + DenseMatrix dshape1_ps, dshape2_ps; + Vector nor; + Vector nL1, nL2; + Vector nM1, nM2; + Vector dshape1_dnM, dshape2_dnM; + DenseMatrix jmat; + #endif + + const int dim = el1.GetDim(); + const int ndofs1 = el1.GetDof(); + const int ndofs2 = (Trans.Elem2No >= 0) ? el2.GetDof() : 0; + const int nvdofs = dim*(ndofs1 + ndofs2); + + // Initially 'elmat' corresponds to the term: + // < { sigma(u) . n }, [v] > = + // < { (lambda div(u) I + mu (grad(u) + grad(u)^T)) . n }, [v] > + // But eventually, it's going to be replaced by: + // elmat := -elmat + alpha*elmat^T + jmat + elmat.SetSize(nvdofs); + elmat = 0.; + + const bool kappa_is_nonzero = (kappa != 0.0); + if (kappa_is_nonzero) + { + jmat.SetSize(nvdofs); + jmat = 0.; + } + + adjJ.SetSize(dim); + shape1.SetSize(ndofs1); + dshape1.SetSize(ndofs1, dim); + dshape1_ps.SetSize(ndofs1, dim); + nor.SetSize(dim); + nL1.SetSize(dim); + nM1.SetSize(dim); + dshape1_dnM.SetSize(ndofs1); + + if (ndofs2) + { + shape2.SetSize(ndofs2); + dshape2.SetSize(ndofs2, dim); + dshape2_ps.SetSize(ndofs2, dim); + nL2.SetSize(dim); + nM2.SetSize(dim); + dshape2_dnM.SetSize(ndofs2); + } + + const IntegrationRule *ir = IntRule; + if (ir == NULL) + { + // a simple choice for the integration order; is this OK? + const int order = 2 * max(el1.GetOrder(), ndofs2 ? el2.GetOrder() : 0); + ir = &IntRules.Get(Trans.GetGeometryType(), order); + } + + //for (int pind = 0; pind < 1; ++pind) + for (int pind = 0; pind < ir->GetNPoints(); ++pind) + { + const IntegrationPoint &ip = ir->IntPoint(pind); + + // Set the integration point in the face and the neighboring elements + Trans.SetAllIntPoints(&ip); + + // Access the neighboring elements' integration points + // Note: eip2 will only contain valid data if Elem2 exists + const IntegrationPoint &eip1 = Trans.GetElement1IntPoint(); + const IntegrationPoint &eip2 = Trans.GetElement2IntPoint(); + + el1.CalcShape(eip1, shape1); + el1.CalcDShape(eip1, dshape1); + + CalcAdjugate(Trans.Elem1->Jacobian(), adjJ); + Mult(dshape1, adjJ, dshape1_ps); + + if (dim == 1) + { + nor(0) = 2*eip1.x - 1.0; + } + else + { + CalcOrtho(Trans.Jacobian(), nor); + } + + double w, wLM; + if (ndofs2) + { + el2.CalcShape(eip2, shape2); + el2.CalcDShape(eip2, dshape2); + CalcAdjugate(Trans.Elem2->Jacobian(), adjJ); + Mult(dshape2, adjJ, dshape2_ps); + + w = ip.weight/2; + const double w2 = w / Trans.Elem2->Weight(); + const double wL2 = w2 * lambda->Eval(*Trans.Elem2, eip2); + const double wM2 = w2 * mu->Eval(*Trans.Elem2, eip2); + nL2.Set(wL2, nor); + nM2.Set(wM2, nor); + wLM = (wL2 + 2.0*wM2); + dshape2_ps.Mult(nM2, dshape2_dnM); + } + else + { + w = ip.weight; + wLM = 0.0; + } + + { + const double w1 = w / Trans.Elem1->Weight(); + const double wL1 = w1 * lambda->Eval(*Trans.Elem1, eip1); + const double wM1 = w1 * mu->Eval(*Trans.Elem1, eip1); + nL1.Set(wL1, nor); + nM1.Set(wM1, nor); + wLM += (wL1 + 2.0*wM1); + dshape1_ps.Mult(nM1, dshape1_dnM); + } + + const double jmatcoef = kappa * (nor*nor) * wLM; + cout<<"jmatcoef is: "< dir_bdr(mesh.bdr_attributes.Max()); + dir_bdr = 0; + dir_bdr[0] = 1; // boundary attribute 1 is Dirichlet + dir_bdr[1] = 1; // boundary attribute 2 is Dirichlet + + BilinearForm a1(&fespace); + a1.AddDomainIntegrator(new ElasticityIntegrator(lambda_c, mu_c)); + a1.AddInteriorFaceIntegrator( + new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa)); + a1.AddBdrFaceIntegrator( + new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr); + a1.Assemble(); + cout<<"a1.Height() is: "< unif(lower_bound, + upper_bound); + default_random_engine re; + + for (size_t i = 0; i < x.Size(); i++) + { + x[i] = unif(re); + //x[i] = 1.0; + //x[i] = 0.0; + //x[i] = i; + } + + y1.SetSize(fespace.GetTrueVSize()); + y1 = 0.0; + + y2.SetSize(fespace.GetTrueVSize()); + y2 = 0.0; + + a1.Mult(x, y1); + a2.Mult(x, y2); + + double norm_diff = 0.0; + cout << "Linear residual norm: " << y1.Norml2() << endl; + cout << "Nonlinear residual norm: " << y2.Norml2() << endl; + + y1 -= y2; + norm_diff = y1.Norml2(); + cout << "Residual difference norm: " << norm_diff << endl; + + a1.Mult(x, y1); + a2.Mult(x, y2); + + y1/=y1.Norml2(); + y2/=y2.Norml2(); + + cout << "Scaled Linear residual norm: " << y1.Norml2() << endl; + cout << "Scaled Nonlinear residual norm: " << y2.Norml2() << endl; + + y1 -= y2; + norm_diff = y1.Norml2(); + cout << "Scaled Residual difference norm: " << norm_diff << endl; + + Operator *J_op = &(a2.GetGradient(x)); + SparseMatrix *J = dynamic_cast(J_op); + + SparseMatrix diff_matrix(*J); + _PrintMatrix(*(diff_matrix.ToDenseMatrix()), "test_elmat1.txt"); + SparseMatrix a1view(a1.SpMat()); + a1view.Finalize(); + _PrintMatrix(*(a1view.ToDenseMatrix()), "test_elmat2.txt"); + + diff_matrix.Add(-1.0, a1.SpMat()); + + norm_diff = diff_matrix.MaxNorm(); + + cout << "Nonlinear Stiffness matrix norm: " << J->MaxNorm() << endl; + cout << "Linear Stiffness matrix norm: " << a1.SpMat().MaxNorm() << endl; + cout << "Stiffness matrix difference norm: " << norm_diff << endl; + + LinearForm b1(&fespace); + b1.AddBdrFaceIntegrator( + new DGElasticityDirichletLFIntegrator( + init_x, lambda_c, mu_c, alpha, kappa), dir_bdr); + b1.Assemble(); + + LinearForm b2(&fespace); + b2.AddBdrFaceIntegrator( + new DGHyperelasticDirichletLFIntegrator( + init_x, &model, alpha, kappa), dir_bdr); + b2.Assemble(); + + cout << "Linear RHS norm: " << b1.Norml2() << endl; + cout << "Nonlinear RHS norm: " << b2.Norml2() << endl; + + b1 -= b2; + norm_diff = b1.Norml2(); + + // Print the norm of the difference + cout << "RHS difference norm: " << norm_diff << endl; + + return; +} + + void ExactSolutionNeoHooke(const Vector &X, Vector &U) + { + int dim = 2; + U = 0.0; + for (size_t i = 0; i < U.Size()/dim; i++) + { + U(i*dim) = pow(X(i*dim), 2.0) + X(i*dim); + U(i*dim+1) = pow(X(i*dim + 1), 2.0) + X(i*dim + 1); + } + } + + void ExactRHSNeoHooke(const Vector &X, Vector &U) + { + int dim = 2; + double K = 3.14; + double mu = 2.33; + U = 0.0; + for (size_t i = 0; i < U.Size()/dim; i++) + { + const double x_1 = X(i*dim); + const double x_2 = X(i*dim + 1); + U(i*dim) = (128.0*K + 128.0*mu + 1024.0*K*x_1 + 1024.0*K*x_2 + 128.0*mu*x_1 + 384.0*mu*x_2 + 3072.0*K*(pow(x_1,2)) + 8192.0*K*x_1*x_2 + 3072.0*K*(pow(x_2,2)) + 128.0*mu*(pow(x_1,2)) + 384.0*mu*(pow(x_2,2)) + 4096.0*K*(pow(x_1,3)) + 24576.0*K*(pow(x_1,2))*x_2 + 24576.0*K*x_1*(pow(x_2,2)) + 4096.0*K*(pow(x_2,3)) + 2048.0*K*(pow(x_1,4)) + 32768.0*K*(pow(x_1,3))*x_2 + 73728.0*K*(pow(x_1,2))*(pow(x_2,2)) + 32768.0*K*x_1*(pow(x_2,3)) + 2048.0*K*(pow(x_2,4)) + 16384.0*K*(pow(x_1,4))*x_2 + 98304.0*K*(pow(x_1,3))*(pow(x_2,2)) + 98304.0*K*(pow(x_1,2))*(pow(x_2,3)) + 16384.0*K*x_1*(pow(x_2,4)) + 49152.0*K*(pow(x_1,4))*(pow(x_2,2)) + 131072.0*K*(pow(x_1,3))*(pow(x_2,3)) + 49152.0*K*(pow(x_1,2))*(pow(x_2,4)) + 65536.0*K*(pow(x_1,4))*(pow(x_2,3)) + 65536.0*K*(pow(x_1,3))*(pow(x_2,4)) + 32768.0*K*(pow(x_1,4))*(pow(x_2,4))) / (64.0 * (pow((1 + 2*x_1), 4))*(pow((1 + 2*x_2),2))); + U(i*dim+1) = (128.0*K + 128.0*mu + 1024.0*K*x_1 + 1024.0*K*x_2 + 384.0*mu*x_1 + 128.0*mu*x_2 + 3072.0*K*(pow(x_1,2)) + 8192.0*K*x_1*x_2 + 3072.0*K*(pow(x_2,2)) + 384.0*mu*(pow(x_1,2)) + 128.0*mu*(pow(x_2,2)) + 4096.0*K*(pow(x_1,3)) + 24576.0*K*(pow(x_1,2))*x_2 + 24576.0*K*x_1*(pow(x_2,2)) + 4096.0*K*(pow(x_2,3)) + 2048.0*K*(pow(x_1,4)) + 32768.0*K*(pow(x_1,3))*x_2 + 73728.0*K*(pow(x_1,2))*(pow(x_2,2)) + 32768.0*K*x_1*(pow(x_2,3)) + 2048.0*K*(pow(x_2,4)) + 16384.0*K*(pow(x_1,4))*x_2 + 98304.0*K*(pow(x_1,3))*(pow(x_2,2)) + 98304.0*K*(pow(x_1,2))*(pow(x_2,3)) + 16384.0*K*x_1*(pow(x_2,4)) + 49152.0*K*(pow(x_1,4))*(pow(x_2,2)) + 131072.0*K*(pow(x_1,3))*(pow(x_2,3)) + 49152.0*K*(pow(x_1,2))*(pow(x_2,4)) + 65536.0*K*(pow(x_1,4))*(pow(x_2,3)) + 65536.0*K*(pow(x_1,3))*(pow(x_2,4)) + 32768.0*K*(pow(x_1,4))*(pow(x_2,4))) / (64.0 * (pow((1 + 2*x_1),2))*(pow((1 + 2*x_2),4))); + } + + //u *= -1.0; + } + + void SimpleExactSolutionNeoHooke(const Vector &X, Vector &U) + { + int dim = 2; + int dof = U.Size()/dim; + U = 0.0; + for (size_t i = 0; i < U.Size()/dim; i++) + { + U(i) = pow(X( i), 2.0) + X( i); + U(dof + i) = pow(X(dof + i), 2.0) + X(dof + i); + } + } + + void SimpleExactRHSNeoHooke(const Vector &X, Vector &U) + { + int dim = 2; + int dof = U.Size()/dim; + + double K = 3.14; + //double mu = 2.33; +/* K = 2.33; + mu = 3.14; */ + U = 0.0; + for (size_t i = 0; i < U.Size()/dim; i++) + { + U(i) = K; + U(dof + i) = K/2.0; + } + + + //u *= -1.0; + } + +// Currently Domain test and Boundary test works. Todo: RHS, Boundary gradients +TEST(TempNeoHookeanStiffnessMatrices, Test_NLElast) +{ + // Test that the nonlinear operators do the correct things + Mesh mesh("meshes/test.2x1.mesh", 1, 1); + + int dim = mesh.Dimension(); + int order = 1; + double alpha = 0.0; // IIPG + double kappa = -1.0; + DG_FECollection fec(order, dim, BasisType::GaussLobatto); + FiniteElementSpace fespace(&mesh, &fec, dim); + + VectorFunctionCoefficient init_x(dim, InitDisplacement); + + Vector lambda(mesh.attributes.Max()); + double _lambda = 3.14; + lambda = _lambda; // Set lambda for all element attributes. + PWConstCoefficient lambda_c(lambda); + Vector mu(mesh.attributes.Max()); + double _mu = 2.33; + //_mu = 0.0; + mu = _mu; // Set mu = 1 for all element attributes. + PWConstCoefficient mu_c(mu); + + Array dir_bdr(mesh.bdr_attributes.Max()); + dir_bdr = 0; + dir_bdr[0] = 1; // boundary attribute 1 is Dirichlet + dir_bdr[1] = 1; // boundary attribute 2 is Dirichlet + + NeoHookeanHypModel model1(_mu, _lambda); + NonlinearForm a1(&fespace); + a1.AddDomainIntegrator(new HyperelasticNLFIntegratorHR(&model1)); + + NeoHookeanModel model2(_mu, _lambda); + NonlinearForm a2(&fespace); + a2.AddDomainIntegrator(new HyperelasticNLFIntegrator(&model2)); + + // Create vectors to hold the values of the forms + Vector x, y0, y1, y2, y3; + + GridFunction x_ref(&fespace); + mesh.GetNodes(x_ref); + x.SetSize(fespace.GetTrueVSize()); + x = x_ref.GetTrueVector(); + + double lower_bound = -1; + double upper_bound = 1; + + //uniform_real_distribution unif(lower_bound, + // upper_bound); + //default_random_engine re; + /* cout<<"couting xi"<(J_op1); + + Operator *J_op2 = &(a2.GetGradient(x)); + SparseMatrix *J2 = dynamic_cast(J_op2); + + //_PrintMatrix(*(diff_matrix.ToDenseMatrix()), "test_elmat1.txt"); + //_PrintMatrix(*(a1view.ToDenseMatrix()), "test_elmat2.txt"); + + cout << "Scaleup NeoHooke matrix norm: " << J1->MaxNorm() << endl; + cout << "MFEM NeoHooke matrix norm: " << J2->MaxNorm() << endl; + + J1->Add(-1.0, *J2); + cout << "Stiffness matrix difference norm: " << J1->MaxNorm() << endl; + + return; +} + + +int main(int argc, char *argv[]) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} \ No newline at end of file diff --git a/test/nonlinear_integ_grad.cpp b/test/nonlinear_integ_grad.cpp index 660c26b6..ad17bbd1 100644 --- a/test/nonlinear_integ_grad.cpp +++ b/test/nonlinear_integ_grad.cpp @@ -5,6 +5,7 @@ #include #include "nonlinear_integ.hpp" #include "hyperreduction_integ.hpp" +#include "nlelast_integ.hpp" #include "input_parser.hpp" #include "etc.hpp" @@ -121,7 +122,8 @@ void CheckGradient(NonlinearFormIntegrator *integ, const IntegratorType type, bo printf("%10s\t%10s\t%10s\t%10s\t%10s\n", "amp", "J0", "J1", "dJdx", "error"); for (int k = 0; k < 40; k++) { - double amp = pow(10.0, -0.25 * k); + //double amp = pow(10.0, -0.25 * k); + double amp = pow(10.0, -5.0-0.25 * k); double dx = amp; if (gg > 1.0e-14) dx /= sqrt(gg); @@ -258,6 +260,57 @@ TEST(TemamTrilinearFormIntegrator, Test_grad) return; } +TEST(HyperelasticNLFIntegratorHR, Test_grad) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["discretization"]["order"] = 1; + + double mu = 2.33; + double K = 3.14; + NeoHookeanHypModel model(mu, K); + //LinElastMaterialModel model(mu, K); + + auto *nlc_nlfi = new HyperelasticNLFIntegratorHR(&model); + + CheckGradient(nlc_nlfi, IntegratorType::DOMAIN, false); + + return; +} + +TEST(DGHyperelasticNLFIntegrator, Test_grad_bdr) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["discretization"]["order"] = 1; + + double mu = 2.33; + double K = 3.14; + NeoHookeanHypModel model(mu, K); + //LinElastMaterialModel model(mu, K); + + auto *nlc_nlfi = new DGHyperelasticNLFIntegrator(&model, 0.0, -1.0); + + CheckGradient(nlc_nlfi, IntegratorType::BDR, true); + + return; +} + +TEST(DGHyperelasticNLFIntegrator, Test_grad_int) +{ + config = InputParser("inputs/dd_mms.yml"); + config.dict_["discretization"]["order"] = 1; + + double mu = 2.33; + double K = 3.14; + NeoHookeanHypModel model(mu, K); + //LinElastMaterialModel model(mu, K); + + auto *nlc_nlfi = new DGHyperelasticNLFIntegrator(&model, 0.0, -1.0); + + CheckGradient(nlc_nlfi, IntegratorType::INTERIOR, true); + + return; +} + int main(int argc, char* argv[]) { MPI_Init(&argc, &argv);