diff --git a/fem/fespace.cpp b/fem/fespace.cpp index c06fb1cf67..2f8e874cd7 100644 --- a/fem/fespace.cpp +++ b/fem/fespace.cpp @@ -1490,10 +1490,8 @@ int FiniteElementSpace::GetNConformingDofs() const int FiniteElementSpace::GetVectorDim() const { const FiniteElement *fe = GetTypicalFE(); - if (fe == nullptr) - { - fe = GetTypicalTraceElement(); - } + MFEM_VERIFY(fe, "A typical finite element does not exist!"); + if (fe->GetRangeType() == FiniteElement::SCALAR) { return GetVDim(); @@ -1501,6 +1499,31 @@ int FiniteElementSpace::GetVectorDim() const return GetVDim()*std::max(GetMesh()->SpaceDimension(), fe->GetRangeDim()); } +int FiniteElementSpace::GetBdrVectorDim() const +{ + const FiniteElement *be = GetTypicalBE(); + MFEM_VERIFY(be, "A typical boundary finite element does not exist!"); + + if (be->GetRangeType() == FiniteElement::SCALAR) + { + return GetVDim(); + } + return GetVDim()*std::max(GetMesh()->SpaceDimension()-1, be->GetRangeDim()); +} + +int FiniteElementSpace::GetFaceVectorDim() const +{ + const FiniteElement *face_el = GetTypicalFaceElement(); + MFEM_VERIFY(face_el, "A typical face finite element does not exist!"); + + if (face_el->GetRangeType() == FiniteElement::SCALAR) + { + return GetVDim(); + } + return GetVDim()*std::max(GetMesh()->SpaceDimension()-1, + face_el->GetRangeDim()); +} + int FiniteElementSpace::GetCurlDim() const { const FiniteElement *fe = GetTypicalFE(); @@ -3915,6 +3938,16 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const return BE; } +const FiniteElement *FiniteElementSpace::GetTypicalBE() const +{ + if (mesh->GetNBE() > 0) { return GetBE(0); } + + Geometry::Type geom = mesh->GetTypicalFaceGeometry(); + const FiniteElement *be = fec->FiniteElementForGeometry(geom); + MFEM_VERIFY(be != nullptr, "Could not determine a typical BE!"); + return be; +} + const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const { MFEM_VERIFY(!IsVariableOrder(), "not implemented"); @@ -3945,6 +3978,11 @@ const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const return fe; } +const FiniteElement *FiniteElementSpace::GetTypicalFaceElement() const +{ + return fec->FiniteElementForGeometry(mesh->GetTypicalFaceGeometry()); +} + const FiniteElement *FiniteElementSpace::GetEdgeElement(int i, int variant) const { diff --git a/fem/fespace.hpp b/fem/fespace.hpp index 3d8fc680d2..1fd395386e 100644 --- a/fem/fespace.hpp +++ b/fem/fespace.hpp @@ -868,6 +868,14 @@ public: */ int GetVectorDim() const; + /// Return the total dimension of a vector on boundary + /** @see GetVectorDim() */ + int GetBdrVectorDim() const; + + /// Return the total dimension of a vector on a face + /** @see GetVectorDim() */ + int GetFaceVectorDim() const; + /// Return the dimension of the curl of a GridFunction defined on this space. /** Note: This assumes a space dimension of 2 or 3 only. */ int GetCurlDim() const; @@ -1333,12 +1341,24 @@ public: associated with i'th boundary face in the mesh object. */ const FiniteElement *GetBE(int i) const; + /// @brief Return a "typical" boundary element. + /// + /// This can be used in situations where the local mesh partition may be + /// empty. + const FiniteElement *GetTypicalBE() const; + /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the mesh object. Faces in this case refer to the MESHDIM-1 primitive so in 2D they are segments and in 1D they are points.*/ const FiniteElement *GetFaceElement(int i) const; + /// @brief Return a "typical" face element. + /// + /// This can be used in situations where the local mesh partition may be + /// empty. + const FiniteElement *GetTypicalFaceElement() const; + /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the mesh object. */ const FiniteElement *GetEdgeElement(int i, int variant = 0) const; diff --git a/fem/gridfunc.cpp b/fem/gridfunc.cpp index 06f700aa84..0cd704dc1b 100644 --- a/fem/gridfunc.cpp +++ b/fem/gridfunc.cpp @@ -351,16 +351,6 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi, } } -int GridFunction::VectorDim() const -{ - return fes->GetVectorDim(); -} - -int GridFunction::CurlDim() const -{ - return fes->GetCurlDim(); -} - void GridFunction::GetTrueDofs(Vector &tv) const { const SparseMatrix *R = fes->GetRestrictionMatrix(); @@ -2660,7 +2650,8 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff, void GridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff, const Array &attr) { - MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()"); + MFEM_VERIFY(BdrVectorDim() == vcoeff.GetVDim(), + "vcoeff vdim != BdrVectorDim()"); Array values_counter; AccumulateAndCountBdrValues(NULL, &vcoeff, attr, values_counter); ComputeMeans(ARITHMETIC, values_counter); @@ -2709,9 +2700,9 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], void GridFunction::ProjectBdrCoefficientNormal( VectorCoefficient &vcoeff, const Array &bdr_attr) { - MFEM_VERIFY(fes->GetVDim() * fes->GetMesh()->SpaceDimension() == - vcoeff.GetVDim(), - "vcoeff vdim mismatch"); + MFEM_VERIFY(BdrVectorDim() == 1, "BdrVectorDim() != 1"); + MFEM_VERIFY(vcoeff.GetVDim() == fes->GetMesh()->SpaceDimension(), + "vcoeff vdim != space dim"); #if 0 // implementation for the case when the face dofs are integrals of the // normal component. @@ -2787,7 +2778,8 @@ void GridFunction::ProjectBdrCoefficientNormal( void GridFunction::ProjectBdrCoefficientTangent( VectorCoefficient &vcoeff, const Array &bdr_attr) { - MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()"); + MFEM_VERIFY(BdrVectorDim()+1 == vcoeff.GetVDim(), + "vcoeff vdim != BdrVectorDim()+1"); Array values_counter; AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter); ComputeMeans(ARITHMETIC, values_counter); diff --git a/fem/gridfunc.hpp b/fem/gridfunc.hpp index d05a787846..d7a0e58ba5 100644 --- a/fem/gridfunc.hpp +++ b/fem/gridfunc.hpp @@ -123,11 +123,21 @@ public: FiniteElementCollection *OwnFEC() { return fec_owned; } - /// Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying #fes - int VectorDim() const; + /** @brief Shortcut for calling FiniteElementSpace::GetVectorDim() on the + underlying #fes */ + int VectorDim() const { return fes->GetVectorDim(); } - /// Shortcut for calling FiniteElementSpace::GetCurlDim() on the underlying #fes - int CurlDim() const; + /** @brief Shortcut for calling FiniteElementSpace::GetBdrVectorDim() on the + underlying #fes */ + int BdrVectorDim() const { return fes->GetBdrVectorDim(); } + + /** @brief Shortcut for calling FiniteElementSpace::GetFaceVectorDim() on + the underlying #fes */ + int FaceVectorDim() const { return fes->GetFaceVectorDim(); } + + /** @brief Shortcut for calling FiniteElementSpace::GetCurlDim() on the + underlying #fes */ + int CurlDim() const { return fes->GetCurlDim(); } /// Read only access to the (optional) internal true-dof Vector. const Vector &GetTrueVector() const diff --git a/fem/pgridfunc.cpp b/fem/pgridfunc.cpp index 3937ff779c..f12aec2eb8 100644 --- a/fem/pgridfunc.cpp +++ b/fem/pgridfunc.cpp @@ -740,14 +740,16 @@ void ParGridFunction::ProjectBdrCoefficient( void ParGridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff, const Array &attr) { - MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()"); + MFEM_VERIFY(BdrVectorDim() == vcoeff.GetVDim(), + "vcoeff vdim != BdrVectorDim()"); ProjectBdrCoefficient(NULL, &vcoeff, attr); } void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, const Array &bdr_attr) { - MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()"); + MFEM_VERIFY(BdrVectorDim()+1 == vcoeff.GetVDim(), + "vcoeff vdim != BdrVectorDim()+1"); Array values_counter; AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter); diff --git a/tests/unit/fem/test_var_order.cpp b/tests/unit/fem/test_var_order.cpp index 50f9616f34..e789b0410e 100644 --- a/tests/unit/fem/test_var_order.cpp +++ b/tests/unit/fem/test_var_order.cpp @@ -696,7 +696,14 @@ static void TestSolveVec(FiniteElementSpace &fespace) GridFunction x(&fespace); x = 0.0; - x.ProjectBdrCoefficient(exsol, ess_attr); + if (x.BdrVectorDim() == 1) + { + x.ProjectBdrCoefficientNormal(exsol, ess_attr); + } + else + { + x.ProjectBdrCoefficientTangent(exsol, ess_attr); + } // Assemble the linear form LinearForm lf(&fespace); @@ -1081,7 +1088,14 @@ static void TestSolveParVec(ParFiniteElementSpace &fespace) ParGridFunction x(&fespace); x = 0.0; - x.ProjectBdrCoefficient(exsol, ess_attr); + if (x.BdrVectorDim() == 1) + { + x.ProjectBdrCoefficientNormal(exsol, ess_attr); + } + else + { + x.ProjectBdrCoefficientTangent(exsol, ess_attr); + } // Assemble the linear form ParLinearForm lf(&fespace); diff --git a/tests/unit/mesh/test_psubmesh.cpp b/tests/unit/mesh/test_psubmesh.cpp index 85ffe94bb5..6efc56a638 100644 --- a/tests/unit/mesh/test_psubmesh.cpp +++ b/tests/unit/mesh/test_psubmesh.cpp @@ -486,8 +486,14 @@ void multidomain_test_3d(FECType fec_type) { cylinder_gf.ProjectCoefficient(vcoeff); outer_gf.ProjectCoefficient(vcoeff); - outer_gf.ProjectBdrCoefficient(vzerocoeff, - outer_cyl_surf_marker); + if (fec_type == FECType::RT) + { + outer_gf.ProjectBdrCoefficientNormal(vzerocoeff, outer_cyl_surf_marker); + } + else + { + outer_gf.ProjectBdrCoefficientTangent(vzerocoeff, outer_cyl_surf_marker); + } outer_gf_ex.ProjectCoefficient(vcoeff); } ParSubMesh::Transfer(cylinder_gf, outer_gf); @@ -507,8 +513,14 @@ void multidomain_test_3d(FECType fec_type) { outer_gf.ProjectCoefficient(vcoeff); cylinder_gf.ProjectCoefficient(vcoeff); - cylinder_gf.ProjectBdrCoefficient(vzerocoeff, - cylinder_cyl_surf_marker); + if (fec_type == FECType::RT) + { + cylinder_gf.ProjectBdrCoefficientNormal(vzerocoeff, cylinder_cyl_surf_marker); + } + else + { + cylinder_gf.ProjectBdrCoefficientTangent(vzerocoeff, cylinder_cyl_surf_marker); + } cylinder_gf_ex.ProjectCoefficient(vcoeff); } ParSubMesh::Transfer(outer_gf, cylinder_gf);