diff --git a/CMakeLists.txt b/CMakeLists.txt index 45745c9f..bffe840f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -200,6 +200,9 @@ set(scaleupROMObj_SOURCES include/poisson_solver.hpp src/poisson_solver.cpp + include/linelast_solver.hpp + src/linelast_solver.cpp + include/stokes_solver.hpp src/stokes_solver.cpp diff --git a/include/interfaceinteg.hpp b/include/interfaceinteg.hpp index 135cbb9a..43921e2f 100644 --- a/include/interfaceinteg.hpp +++ b/include/interfaceinteg.hpp @@ -11,212 +11,244 @@ namespace mfem { -class InterfaceNonlinearFormIntegrator : virtual public HyperReductionIntegrator -{ -protected: - InterfaceNonlinearFormIntegrator(const bool precomputable_ = false, const IntegrationRule *ir = NULL) - : HyperReductionIntegrator(precomputable_, ir) {} -public: - // FaceElementTransformations belongs to one mesh (having mesh pointer). - // In order to extract element/transformation from each mesh, - // two FaceElementTransformations are needed. - virtual void AssembleInterfaceVector(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Tr1, - FaceElementTransformations &Tr2, - const Vector &elfun1, const Vector &elfun2, - Vector &elvect1, Vector &elvect2); - - /// @brief Assemble the local action of the gradient of the - /// NonlinearFormIntegrator resulting from a face integral term. - virtual void AssembleInterfaceGrad(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Tr1, - FaceElementTransformations &Tr2, - const Vector &elfun1, const Vector &elfun2, - Array2D &elmats); - - virtual void AssembleInterfaceMatrix(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Trans1, - FaceElementTransformations &Trans2, - Array2D &elmats); - - /** Abstract method used for assembling InteriorFaceIntegrators in a - MixedBilinearFormDGExtension. */ - virtual void AssembleInterfaceMatrix(const FiniteElement &trial_fe1, - const FiniteElement &trial_fe2, - const FiniteElement &test_fe1, - const FiniteElement &test_fe2, - FaceElementTransformations &Trans1, - FaceElementTransformations &Trans2, - Array2D &elmats) - { mfem_error("Abstract method InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix!\n"); } -}; - -class InterfaceDGDiffusionIntegrator : public InterfaceNonlinearFormIntegrator -{ -protected: - Coefficient *Q; - MatrixCoefficient *MQ; - double sigma, kappa; - - // these are not thread-safe! - Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni; - DenseMatrix dshape1, dshape2, mq, adjJ; - Vector nor2; - Array2D jmats; - -public: - InterfaceDGDiffusionIntegrator(const double s, const double k) - : Q(NULL), MQ(NULL), sigma(s), kappa(k) { } - InterfaceDGDiffusionIntegrator(Coefficient &q, const double s, const double k) - : Q(&q), MQ(NULL), sigma(s), kappa(k) { } - InterfaceDGDiffusionIntegrator(MatrixCoefficient &q, const double s, const double k) - : Q(NULL), MQ(&q), sigma(s), kappa(k) { } - - virtual void AssembleInterfaceMatrix(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Trans1, - FaceElementTransformations &Trans2, - Array2D &elmats); -}; - -// DGDiffusionFaceIntegrator -class InterfaceDGVectorDiffusionIntegrator : public InterfaceNonlinearFormIntegrator -{ -public: - InterfaceDGVectorDiffusionIntegrator(double alpha_, double kappa_) - : mu(NULL), alpha(alpha_), kappa(kappa_) { } - - InterfaceDGVectorDiffusionIntegrator(Coefficient &mu_, - double alpha_, double kappa_) - : mu(&mu_), alpha(alpha_), kappa(kappa_) { } - - virtual void AssembleInterfaceMatrix(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Trans1, - FaceElementTransformations &Trans2, - Array2D &elmats); - - using InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix; - -protected: - Coefficient *mu; - double alpha, kappa; - - // 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 nM1, nM2; // nM1 = (mu1 * ip.weight / detJ1) nor - Vector dshape1_dnM, dshape2_dnM; // dshape1_dnM = dshape1_ps . nM1 - // 'jmat' corresponds to the term: kappa - Array2D jmats; - - // Since elmats are already blocked out, we do not need the offsets. - // offsets are used only for jmat, to determine the lower-triangular part. - 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 &row_shape, const Vector &col_shape, const Vector &col_dshape_dnM, - DenseMatrix &elmat, DenseMatrix &jmat); -}; - -class InterfaceDGNormalFluxIntegrator : public InterfaceNonlinearFormIntegrator -{ -private: - int dim; - int order; - int p; - - int trial_dof1, trial_dof2, test_dof1, test_dof2; - int trial_vdof1, trial_vdof2; - - double w, wn; - int i, j, idof, jdof, jm; - - Vector nor, wnor; - Vector shape1, shape2; - // Vector divshape; - Vector trshape1, trshape2; - // DenseMatrix vshape1, vshape2; - // Vector vshape1_n, vshape2_n; - -public: - InterfaceDGNormalFluxIntegrator() {}; - virtual ~InterfaceDGNormalFluxIntegrator() {}; - - virtual void AssembleInterfaceMatrix(const FiniteElement &trial_fe1, - const FiniteElement &trial_fe2, - const FiniteElement &test_fe1, - const FiniteElement &test_fe2, - FaceElementTransformations &Trans1, - FaceElementTransformations &Trans2, - Array2D &elmats); - - using InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix; -}; - -class InterfaceDGTemamFluxIntegrator : public InterfaceNonlinearFormIntegrator -{ -private: - int dim, ndofs1, ndofs2, nvdofs1, nvdofs2; - double w, un; - Coefficient *Q{}; - - Vector nor, shape1, shape2, u1, u2, flux; - DenseMatrix udof1, udof2, elv1, elv2; - DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21; - -public: - InterfaceDGTemamFluxIntegrator(Coefficient &q) : Q(&q) {}; - virtual ~InterfaceDGTemamFluxIntegrator() {}; - - virtual void AssembleInterfaceVector(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Tr1, - FaceElementTransformations &Tr2, - const Vector &elfun1, const Vector &elfun2, - Vector &elvect1, Vector &elvect2); - - virtual void AssembleInterfaceGrad(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Tr1, - FaceElementTransformations &Tr2, - const Vector &elfun1, const Vector &elfun2, - Array2D &elmats); -}; - -// DGElasticityIntegrator -// TODO(axel): you'll need to implement this class for dg interface for two different meshes. -/* This is pretty much the copy version of MFEM DGElasticityIntegrator class, - which works for two different meshes. - For reference, the Poisson equation uses InterfaceDGDiffusionIntegrator, - which is the copy version of DGDiffusionIntegrator. - */ -class InterfaceDGElasticityIntegrator : public InterfaceNonlinearFormIntegrator -{ -public: - InterfaceDGElasticityIntegrator(){ } - - virtual void AssembleInterfaceMatrix(const FiniteElement &el1, - const FiniteElement &el2, - FaceElementTransformations &Trans1, - FaceElementTransformations &Trans2, - Array2D &elmats) {} - - using InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix; -}; + class InterfaceNonlinearFormIntegrator : virtual public HyperReductionIntegrator + { + protected: + InterfaceNonlinearFormIntegrator(const bool precomputable_ = false, const IntegrationRule *ir = NULL) + : HyperReductionIntegrator(precomputable_, ir) {} + + public: + // FaceElementTransformations belongs to one mesh (having mesh pointer). + // In order to extract element/transformation from each mesh, + // two FaceElementTransformations are needed. + virtual void AssembleInterfaceVector(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Tr1, + FaceElementTransformations &Tr2, + const Vector &elfun1, const Vector &elfun2, + Vector &elvect1, Vector &elvect2); + + /// @brief Assemble the local action of the gradient of the + /// NonlinearFormIntegrator resulting from a face integral term. + virtual void AssembleInterfaceGrad(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Tr1, + FaceElementTransformations &Tr2, + const Vector &elfun1, const Vector &elfun2, + Array2D &elmats); + + virtual void AssembleInterfaceMatrix(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Trans1, + FaceElementTransformations &Trans2, + Array2D &elmats); + + /** Abstract method used for assembling InteriorFaceIntegrators in a + MixedBilinearFormDGExtension. */ + virtual void AssembleInterfaceMatrix(const FiniteElement &trial_fe1, + const FiniteElement &trial_fe2, + const FiniteElement &test_fe1, + const FiniteElement &test_fe2, + FaceElementTransformations &Trans1, + FaceElementTransformations &Trans2, + Array2D &elmats) + { + mfem_error("Abstract method InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix!\n"); + } + }; + + class InterfaceDGDiffusionIntegrator : public InterfaceNonlinearFormIntegrator + { + protected: + Coefficient *Q; + MatrixCoefficient *MQ; + double sigma, kappa; + + // these are not thread-safe! + Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni; + DenseMatrix dshape1, dshape2, mq, adjJ; + Vector nor2; + Array2D jmats; + + public: + InterfaceDGDiffusionIntegrator(const double s, const double k) + : Q(NULL), MQ(NULL), sigma(s), kappa(k) {} + InterfaceDGDiffusionIntegrator(Coefficient &q, const double s, const double k) + : Q(&q), MQ(NULL), sigma(s), kappa(k) {} + InterfaceDGDiffusionIntegrator(MatrixCoefficient &q, const double s, const double k) + : Q(NULL), MQ(&q), sigma(s), kappa(k) {} + + virtual void AssembleInterfaceMatrix(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Trans1, + FaceElementTransformations &Trans2, + Array2D &elmats); + }; + + // DGDiffusionFaceIntegrator + class InterfaceDGVectorDiffusionIntegrator : public InterfaceNonlinearFormIntegrator + { + public: + InterfaceDGVectorDiffusionIntegrator(double alpha_, double kappa_) + : mu(NULL), alpha(alpha_), kappa(kappa_) {} + + InterfaceDGVectorDiffusionIntegrator(Coefficient &mu_, + double alpha_, double kappa_) + : mu(&mu_), alpha(alpha_), kappa(kappa_) {} + + virtual void AssembleInterfaceMatrix(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Trans1, + FaceElementTransformations &Trans2, + Array2D &elmats); + + using InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix; + + protected: + Coefficient *mu; + double alpha, kappa; + + // 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 nM1, nM2; // nM1 = (mu1 * ip.weight / detJ1) nor + Vector dshape1_dnM, dshape2_dnM; // dshape1_dnM = dshape1_ps . nM1 + // 'jmat' corresponds to the term: kappa + Array2D jmats; + + // Since elmats are already blocked out, we do not need the offsets. + // offsets are used only for jmat, to determine the lower-triangular part. + 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 &row_shape, const Vector &col_shape, const Vector &col_dshape_dnM, + DenseMatrix &elmat, DenseMatrix &jmat); + }; + + class InterfaceDGNormalFluxIntegrator : public InterfaceNonlinearFormIntegrator + { + private: + int dim; + int order; + int p; + + int trial_dof1, trial_dof2, test_dof1, test_dof2; + int trial_vdof1, trial_vdof2; + + double w, wn; + int i, j, idof, jdof, jm; + + Vector nor, wnor; + Vector shape1, shape2; + // Vector divshape; + Vector trshape1, trshape2; + // DenseMatrix vshape1, vshape2; + // Vector vshape1_n, vshape2_n; + + public: + InterfaceDGNormalFluxIntegrator(){}; + virtual ~InterfaceDGNormalFluxIntegrator(){}; + + virtual void AssembleInterfaceMatrix(const FiniteElement &trial_fe1, + const FiniteElement &trial_fe2, + const FiniteElement &test_fe1, + const FiniteElement &test_fe2, + FaceElementTransformations &Trans1, + FaceElementTransformations &Trans2, + Array2D &elmats); + + using InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix; + }; + + class InterfaceDGTemamFluxIntegrator : public InterfaceNonlinearFormIntegrator + { + private: + int dim, ndofs1, ndofs2, nvdofs1, nvdofs2; + double w, un; + Coefficient *Q{}; + + Vector nor, shape1, shape2, u1, u2, flux; + DenseMatrix udof1, udof2, elv1, elv2; + DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21; + + public: + InterfaceDGTemamFluxIntegrator(Coefficient &q) : Q(&q) {}; + virtual ~InterfaceDGTemamFluxIntegrator() {}; + + virtual void AssembleInterfaceVector(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Tr1, + FaceElementTransformations &Tr2, + const Vector &elfun1, const Vector &elfun2, + Vector &elvect1, Vector &elvect2); + + virtual void AssembleInterfaceGrad(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Tr1, + FaceElementTransformations &Tr2, + const Vector &elfun1, const Vector &elfun2, + Array2D &elmats); + }; + + // DGElasticityIntegrator + // TODO(axel): you'll need to implement this class for dg interface for two different meshes. + /* This is pretty much the copy version of MFEM DGElasticityIntegrator class, + which works for two different meshes. + For reference, the Poisson equation uses InterfaceDGDiffusionIntegrator, + which is the copy version of DGDiffusionIntegrator. + */ + class InterfaceDGElasticityIntegrator : public InterfaceNonlinearFormIntegrator + { + protected: + double alpha, kappa; + + // 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 nor1, nor2; // nor = |weight(J_face)| n + Vector nM1, nM2; // nM1 = (mu1 * ip.weight / detJ1) nor + Vector nL1, nL2; + Vector dshape1_dnM, dshape2_dnM; // dshape1_dnM = dshape1_ps . nM1 + // 'jmat' corresponds to the term: kappa + Array2D jmats; + + PWConstCoefficient *lambda, *mu; + + public: + InterfaceDGElasticityIntegrator(double alpha_, double kappa_) : alpha(alpha_), kappa(kappa_) {}; + + virtual 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); + + virtual void AssembleInterfaceMatrix(const FiniteElement &el1, + const FiniteElement &el2, + FaceElementTransformations &Trans1, + FaceElementTransformations &Trans2, + Array2D &elmats); + + using InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix; + }; } // namespace mfem diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp index bd86e542..98d4621d 100644 --- a/include/linelast_solver.hpp +++ b/include/linelast_solver.hpp @@ -44,11 +44,11 @@ class LinElastSolver : public MultiBlockSolver Array mu_cs; // DG parameters specific to linear elasticity equation. - double sigma = -1.0; + double alpha = -1.0; double kappa = -1.0; // Initial positions - VectorFunctionCoefficient init_x; + VectorFunctionCoefficient *init_x; public: LinElastSolver(); @@ -68,28 +68,21 @@ class LinElastSolver : public MultiBlockSolver u(u.Size() - 1) = -0.2 * x(0); } - virtual void SetupBCVariables() override; - virtual void AddBCFunction(std::function F, const int battr = -1); - virtual void AddBCFunction(const double &F, const int battr = -1); + //virtual void SetupBCVariables() override; + //virtual void AddBCFunction(std::function F, const int battr = -1); + //virtual void AddBCFunction(const double &F, const int battr = -1); virtual void InitVariables(); + virtual void SetupMaterialVariables(); virtual void BuildOperators(); virtual void BuildRHSOperators(); virtual void BuildDomainOperators(); - virtual bool BCExistsOnBdr(const int &global_battr_idx); - virtual void SetupBCOperators(); - virtual void SetupRHSBCOperators(); - virtual void SetupDomainBCOperators(); + //virtual bool BCExistsOnBdr(const int &global_battr_idx); + //virtual void SetupBCOperators(); + //virtual void SetupRHSBCOperators(); + //virtual void SetupDomainBCOperators(); - virtual void AddRHSFunction(std::function F) - { - rhs_coeffs.Append(new FunctionCoefficient(F)); - } - virtual void AddRHSFunction(const double F) - { - rhs_coeffs.Append(new ConstantCoefficient(F)); - } virtual void Assemble(); virtual void AssembleRHS(); @@ -99,17 +92,17 @@ class LinElastSolver : public MultiBlockSolver virtual void AssembleInterfaceMatrixes(); // Component-wise assembly - virtual void BuildCompROMElement(Array &fes_comp); - virtual void BuildBdrROMElement(Array &fes_comp); - virtual void BuildInterfaceROMElement(Array &fes_comp); + //virtual void BuildCompROMElement(Array &fes_comp); + //virtual void BuildBdrROMElement(Array &fes_comp); + //virtual void BuildInterfaceROMElement(Array &fes_comp); virtual bool Solve(); - virtual void ProjectOperatorOnReducedBasis(); + //virtual void ProjectOperatorOnReducedBasis(); - void SanityCheckOnCoeffs(); + //void SanityCheckOnCoeffs(); - virtual void SetParameterizedProblem(ParameterizedProblem *problem); + //virtual void SetParameterizedProblem(ParameterizedProblem *problem); }; #endif diff --git a/src/interfaceinteg.cpp b/src/interfaceinteg.cpp index 5cda4c7c..57b8b3d5 100644 --- a/src/interfaceinteg.cpp +++ b/src/interfaceinteg.cpp @@ -1014,16 +1014,16 @@ namespace mfem DenseMatrix dshape1, dshape2; DenseMatrix adjJ; DenseMatrix dshape1_ps, dshape2_ps; - Vector nor; + Vector nor1, nor2; Vector nL1, nL2; Vector nM1, nM2; Vector dshape1_dnM, dshape2_dnM; DenseMatrix jmat; #endif - + bool boundary = false; const int dim = el1.GetDim(); const int ndofs1 = el1.GetDof(); - const int ndofs2 = (boundary) ? el2.GetDof() : 0; + const int ndofs2 = (boundary) ? 0 : el2.GetDof(); const int nvdofs = dim * (ndofs1 + ndofs2); // Initially 'elmat' corresponds to the term: @@ -1031,10 +1031,10 @@ namespace mfem // < { (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 - elmats(0, 0)->SetSize(ndof1, ndof1); - elmats(0, 1)->SetSize(ndof1, ndof2); - elmats(1, 0)->SetSize(ndof2, ndof1); - elmats(1, 1)->SetSize(ndof2, ndof2); + elmats(0, 0)->SetSize(ndofs1, ndofs1); + elmats(0, 1)->SetSize(ndofs1, ndofs2); + elmats(1, 0)->SetSize(ndofs2, ndofs1); + elmats(1, 1)->SetSize(ndofs2, ndofs2); for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) *elmats(i, j) = 0.0; @@ -1045,11 +1045,11 @@ namespace mfem if (kappa_is_nonzero) { jmats.SetSize(2, 2); - jmats(0, 0) = new DenseMatrix(ndof1, ndof1); + jmats(0, 0) = new DenseMatrix(ndofs1, ndofs1); // only the lower-triangular part of jmat is assembled. jmats(0, 1) = NULL; - jmats(1, 0) = new DenseMatrix(ndof2, ndof1); - jmats(1, 1) = new DenseMatrix(ndof2, ndof2); + jmats(1, 0) = new DenseMatrix(ndofs2, ndofs1); + jmats(1, 1) = new DenseMatrix(ndofs2, ndofs2); for (int i = 0; i < 2; i++) for (int j = 0; j <= i; j++) *jmats(i, j) = 0.0; @@ -1061,7 +1061,7 @@ namespace mfem shape1.SetSize(ndofs1); dshape1.SetSize(ndofs1, dim); dshape1_ps.SetSize(ndofs1, dim); - nor.SetSize(dim); + nor1.SetSize(dim); nL1.SetSize(dim); nM1.SetSize(dim); dshape1_dnM.SetSize(ndofs1); @@ -1071,6 +1071,7 @@ namespace mfem shape2.SetSize(ndofs2); dshape2.SetSize(ndofs2, dim); dshape2_ps.SetSize(ndofs2, dim); + nor2.SetSize(dim); nL2.SetSize(dim); nM2.SetSize(dim); dshape2_dnM.SetSize(ndofs2); @@ -1081,7 +1082,7 @@ namespace mfem { // 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); + ir = &IntRules.Get(Trans1.GetGeometryType(), order); assert(Trans1.GetGeometryType() == Trans2.GetGeometryType()); } @@ -1101,12 +1102,12 @@ namespace mfem // computing outward normal vectors. if (dim == 1) { - nor(0) = 2 * eip1.x - 1.0; + nor1(0) = 2 * eip1.x - 1.0; nor2(0) = 2 * eip2.x - 1.0; } else { - CalcOrtho(Trans1.Jacobian(), nor); + CalcOrtho(Trans1.Jacobian(), nor1); CalcOrtho(Trans2.Jacobian(), nor2); } @@ -1128,8 +1129,8 @@ namespace mfem const double w2 = w / Trans2.Elem1->Weight(); const double wL2 = w2 * lambda->Eval(*Trans2.Elem1, eip2); const double wM2 = w2 * mu->Eval(*Trans2.Elem1, eip2); - nL2.Set(wL2, nor); - nM2.Set(wM2, nor); + nL2.Set(wL2, nor2); + nM2.Set(wM2, nor2); wLM = (wL2 + 2.0 * wM2); dshape2_ps.Mult(nM2, dshape2_dnM); } @@ -1143,13 +1144,13 @@ namespace mfem const double w1 = w / Trans1.Elem1->Weight(); const double wL1 = w1 * lambda->Eval(*Trans1.Elem1, eip1); const double wM1 = w1 * mu->Eval(*Trans1.Elem1, eip1); - nL1.Set(wL1, nor); - nM1.Set(wM1, nor); + nL1.Set(wL1, nor1); + nM1.Set(wM1, nor1); wLM += (wL1 + 2.0 * wM1); dshape1_ps.Mult(nM1, dshape1_dnM); } - const double jmatcoef = kappa * (nor * nor) * wLM; + const double jmatcoef = kappa * (nor1 * nor2) * wLM; // (1,1) block AssembleBlock( dim, ndofs1, ndofs1, 0, 0, jmatcoef, nL1, nM1, @@ -1179,8 +1180,10 @@ namespace mfem // elmat := -elmat + sigma*elmat^t + jmat Array ndof_array(2); - ndof_array[0] = ndof1; - ndof_array[1] = ndof2; + ndof_array[0] = ndofs1; + ndof_array[1] = ndofs2; + DenseMatrix *elmat = NULL; + DenseMatrix *jmat = NULL; DenseMatrix *elmat12 = NULL; if (kappa_is_nonzero) { @@ -1193,23 +1196,23 @@ namespace mfem for (int j = 0; j < i; j++) { double aij = (*elmat)(i, j), aji = (*elmat)(j, i), mij = (*jmat)(i, j); - (*elmat)(i, j) = sigma * aji - aij + mij; - (*elmat)(j, i) = sigma * aij - aji + mij; + (*elmat)(i, j) = alpha * aji - aij + mij; + (*elmat)(j, i) = alpha * aij - aji + mij; } - (*elmat)(i, i) = (sigma - 1.) * (*elmat)(i, i) + (*jmat)(i, i); + (*elmat)(i, i) = (alpha - 1.) * (*elmat)(i, i) + (*jmat)(i, i); } // for (int i = 0; i < ndofs_array[I]; i++) } // for (int I = 0; I < 2; I++) elmat = elmats(1, 0); jmat = jmats(1, 0); assert(jmat != NULL); elmat12 = elmats(0, 1); - for (int i = 0; i < ndof2; i++) + for (int i = 0; i < ndofs2; i++) { - for (int j = 0; j < ndof1; j++) + for (int j = 0; j < ndofs1; j++) { double aij = (*elmat)(i, j), aji = (*elmat12)(j, i), mij = (*jmat)(i, j); - (*elmat)(i, j) = sigma * aji - aij + mij; - (*elmat12)(j, i) = sigma * aij - aji + mij; + (*elmat)(i, j) = alpha * aji - aij + mij; + (*elmat12)(j, i) = alpha * aij - aji + mij; } // for (int j = 0; j < ndofs1; j++) } // for (int i = 0; i < ndofs2; i++) } // if (kappa_is_nonzero) @@ -1223,21 +1226,21 @@ namespace mfem for (int j = 0; j < i; j++) { double aij = (*elmat)(i, j), aji = (*elmat)(j, i); - (*elmat)(i, j) = sigma * aji - aij; - (*elmat)(j, i) = sigma * aij - aji; + (*elmat)(i, j) = alpha * aji - aij; + (*elmat)(j, i) = alpha * aij - aji; } - (*elmat)(i, i) *= (sigma - 1.); + (*elmat)(i, i) *= (alpha - 1.); } // for (int i = 0; i < ndofs_array[I]; i++) } // for (int I = 0; I < 2; I++) elmat = elmats(1, 0); elmat12 = elmats(0, 1); - for (int i = 0; i < ndof2; i++) + for (int i = 0; i < ndofs2; i++) { - for (int j = 0; j < ndof1; j++) + for (int j = 0; j < ndofs1; j++) { double aij = (*elmat)(i, j), aji = (*elmat12)(j, i); - (*elmat)(i, j) = sigma * aji - aij; - (*elmat12)(j, i) = sigma * aij - aji; + (*elmat)(i, j) = alpha * aji - aij; + (*elmat12)(j, i) = alpha * aij - aji; } // for (int j = 0; j < ndofs1; j++) } // for (int i = 0; i < ndofs2; i++) } // not if (kappa_is_nonzero) @@ -1247,10 +1250,10 @@ namespace mfem } // static method - void InterfaceDGElasticityIntegrator::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 InterfaceDGElasticityIntegrator::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) { // row_offset and col_offset are not needed for elmat. for (int jm = 0, j = 0; jm < dim; ++jm) diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index cd5ce3ba..c9350d69 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -13,7 +13,7 @@ using namespace mfem; LinElastSolver::LinElastSolver() : MultiBlockSolver() { - sigma = config.GetOption("discretization/interface/sigma", -1.0); + alpha = config.GetOption("discretization/interface/alpha", -1.0); kappa = config.GetOption("discretization/interface/kappa", (order + 1) * (order + 1)); var_names = GetVariableNames(); @@ -41,7 +41,7 @@ LinElastSolver::LinElastSolver() fes[m] = new FiniteElementSpace(meshes[m], fec[0], udim); } - init_x = VectorFunctionCoefficient(dim, InitDisplacement); + init_x = new VectorFunctionCoefficient(dim, InitDisplacement); } LinElastSolver::~LinElastSolver() @@ -83,8 +83,8 @@ void LinElastSolver::SetupMaterialVariables() for (int m = 0; m < numSub; m++) { - lambda_cs[m] = lambda_c; - mu_cs[m] = mu_c; + lambda_cs[m] = &lambda_c; + mu_cs[m] = &mu_c; } } @@ -129,7 +129,7 @@ void LinElastSolver::InitVariables() for (int m = 0; m < numSub; m++) { us[m] = new GridFunction(fes[m], U->GetBlock(m), 0); - us[m]->ProjectCoefficient(init_x); + us[m]->ProjectCoefficient(*init_x); // BC's are weakly constrained and there is no essential dofs. // Does this make any difference? @@ -153,7 +153,7 @@ void LinElastSolver::BuildRHSOperators() for (int m = 0; m < numSub; m++) { bs[m] = new LinearForm(fes[m], RHS->GetBlock(m).GetData()); - bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(init_x, lambda_cs[m], mu_cs[m], alpha, kappa), bdr_marker[0]); // TODO bdr_marker will be extended to include forces too + bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*init_x, *(lambda_cs[m]), *(mu_cs[m]), alpha, kappa), *(bdr_markers[0])); //m ight be wrong TODO bdr_marker will be extended to include forces too, also, unsure if this is correct } } @@ -166,18 +166,18 @@ void LinElastSolver::BuildDomainOperators() for (int m = 0; m < numSub; m++) { as[m] = new BilinearForm(fes[m]); - as[m]->AddDomainIntegrator(new ElasticityIntegrator(lambda_cs[m], mu_cs[m])); + as[m]->AddDomainIntegrator(new ElasticityIntegrator(*(lambda_cs[m]), *(mu_cs[m]))); if (full_dg) { as[m]->AddInteriorFaceIntegrator( - new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa)); + new DGElasticityIntegrator(*(lambda_cs[m]), *(mu_cs[m]), alpha, kappa)); as[m]->AddBdrFaceIntegrator( - new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr); + new DGElasticityIntegrator(*(lambda_cs[m]), *(mu_cs[m]), alpha, kappa), *(bdr_markers[0])); } } a_itf = new InterfaceForm(meshes, fes, topol_handler); - a_itf->AddIntefaceIntegrator(new InterfaceDGElasticityIntegrator(sigma, kappa)); + a_itf->AddIntefaceIntegrator(new InterfaceDGElasticityIntegrator(alpha, kappa)); } void LinElastSolver::Assemble() diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 0133db8d..0299ebc6 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -6,6 +6,7 @@ #include "component_topology_handler.hpp" #include "multiblock_solver.hpp" #include "poisson_solver.hpp" +#include "linelast_solver.hpp" #include "stokes_solver.hpp" #include "steady_ns_solver.hpp" #include "etc.hpp" @@ -53,7 +54,7 @@ MultiBlockSolver* InitSolver() if (solver_type == "poisson") { solver = new PoissonSolver; } else if (solver_type == "stokes") { solver = new StokesSolver; } else if (solver_type == "steady-ns") { solver = new SteadyNSSolver; } - else if (solver_type == "linelast") { solver = new LinElastSolver; } + //else if (solver_type == "linelast") { solver = new LinElastSolver; } // TODO: make this work else { printf("Unknown MultiBlockSolver %s!\n", solver_type.c_str());