diff --git a/CMakeLists.txt b/CMakeLists.txt index 45dc805..3030740 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,6 +26,8 @@ set(MESHFIELD_HEADERS src/MeshField_Macros.hpp src/KokkosController.hpp src/MeshField_Element.hpp + src/MeshField_LinearAlgebra.hpp + src/MeshField_MappedElement.hpp src/MeshField_Integrate.hpp src/MeshField_Shape.hpp src/MeshField_ShapeField.hpp @@ -201,40 +203,30 @@ function(meshfields_add_exe EXE_NAME EXE_SRC) endif() endfunction() -meshfields_add_smoke_exe(KokkosTests test/testKokkos.cpp) +meshfields_add_smoke_exe(ControllerTests test/testController.cpp) meshfields_add_exe(SerializationTests test/testSerialize.cpp) -meshfields_add_exe(ElementTests test/testElement.cpp) meshfields_add_exe(ControllerLifetimeTest test/testControllerLifetime.cpp) -meshfields_add_exe(ElementJacobian1d test/testElementJacobian1d.cpp) -meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp) -meshfields_add_exe(ElementJacobian3d test/testElementJacobian3d.cpp) -meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp) -meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp) +meshfields_add_exe(ElementJacobian test/testElementJacobian.cpp) +meshfields_add_exe(Reproduction test/testReproduction.cpp) meshfields_add_exe(ExceptionTest test/testExceptions.cpp) -meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) -meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) +meshfields_add_exe(Integrator test/testIntegrator.cpp) +meshfields_add_exe(Interpolate test/testInterpolate.cpp) +meshfields_add_exe(OmegahMappedElement test/testOmegahMappedElement.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp) - meshfields_add_smoke_exe(CabanaTests test/testCabana.cpp) - meshfields_add_exe(CabanaControllerLifetimeTest test/testCabanaControllerLifetime.cpp) - smoke_test_func(CabanaTests ./CabanaTests) test_func(ControllerPerformance ./ControllerPerformance) - test_func(CabanaControllerLifetimeTest ./CabanaControllerLifetimeTest) endif() -smoke_test_func(KokkosTests ./KokkosTests) +smoke_test_func(ControllerTests ./ControllerTests) test_func(SerializationTests ./SerializationTests) -test_func(ElementTests ./ElementTests) test_func(ControllerLifetimeTest ./ControllerLifetimeTest) -test_func(ElementJacobian1d ./ElementJacobian1d) -test_func(ElementJacobian2d ./ElementJacobian2d) -test_func(ElementJacobian3d ./ElementJacobian3d) -test_func(CountIntegrator ./CountIntegrator) -smoke_test_func(OmegahTriTests ./OmegahTriTests) -test_func(OmegahTetTest ./OmegahTetTest) -test_func(IntegratorMeasure ./IntegratorMeasure) +test_func(ElementJacobian ./ElementJacobian) +test_func(Reproduction ./Reproduction) +test_func(Integrator ./Integrator) +test_func(Interpolate ./Interpolate) +test_func(OmegahMappedElement ./OmegahMappedElement) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index dd541b8..e67421b 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -5,6 +5,7 @@ #include "MeshField_Element.hpp" #include "MeshField_Fail.hpp" #include "MeshField_For.hpp" +#include "MeshField_MappedElement.hpp" #include "MeshField_ShapeField.hpp" #include "Omega_h_file.hpp" //move #include "Omega_h_mesh.hpp" //move @@ -31,30 +32,27 @@ MeshField::MeshInfo getMeshInfo(Omega_h::Mesh &mesh) { } return meshInfo; } - -template - typename Controller = MeshField::KokkosController> -decltype(MeshField::CreateCoordinateField( - MeshField::MeshInfo())) -createCoordinateField(const MeshField::MeshInfo &mesh_info, - Omega_h::Reals coords) { - const auto meshDim = mesh_info.dim; - auto coordFieldWithCtrlr = - MeshField::CreateCoordinateField( - mesh_info); - auto coordField = coordFieldWithCtrlr.field; - auto setCoordField = KOKKOS_LAMBDA(const int &i) { - coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim]; - coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1]; - if constexpr (dim == 3) { - coordField(i, 0, 2, MeshField::Vertex) = coords[i * meshDim + 2]; - } - }; - MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx}, - setCoordField, "setCoordField"); - return coordFieldWithCtrlr; -} +// we use functors because cannot use auto return type with KOKKOS_LAMBDAS +template +struct SetCoordFieldFunctor +{ + SetCoordFieldFunctor(Field& field, Coords& coords, int meshDim, + const CoordinateMap& map) : + field_(field), coords_(coords), meshDim_(meshDim), map_(map) {} + KOKKOS_INLINE_FUNCTION + void operator()(const size_t vtx) const { + Kokkos::Array x{}; + for (size_t d = 0; d < meshDim_; ++d) + x[d] = coords_[vtx * meshDim_ + d]; + const auto mapped = map_(x); + for (size_t d = 0; d < meshDim_; ++d) + field_(vtx, 0, d, MeshField::Vertex) = mapped[d]; + } + Field field_; + Coords coords_; + int meshDim_; + const CoordinateMap map_; +}; } // anonymous namespace @@ -284,6 +282,16 @@ template auto getTetrahedronElement(Omega_h::Mesh &mesh) { } // namespace Omegah +/** + * @brief the CoordinateMap that leaves the mesh geometry unchanged + */ +struct IdentityCoordinateMap { + template + KOKKOS_INLINE_FUNCTION Point operator()(Point const &x) const { + return x; + } +}; + template typename Controller = MeshField::KokkosController> @@ -292,15 +300,14 @@ class OmegahMeshField { Omega_h::Mesh &mesh; const MeshField::MeshInfo meshInfo; using CoordField = - decltype(createCoordinateField( - MeshField::MeshInfo(), Omega_h::Reals())); + decltype(MeshField::CreateCoordinateField(MeshField::MeshInfo())); CoordField coordField; public: OmegahMeshField(Omega_h::Mesh &mesh_in) : mesh(mesh_in), meshInfo(getMeshInfo(mesh)), - coordField(createCoordinateField( - getMeshInfo(mesh_in), mesh_in.coords())) { + coordField(CreateLagrangeCoordinateField<1>()) { static_assert(dim == 1 || dim == 2 || dim == 3); } @@ -311,8 +318,81 @@ class OmegahMeshField { order, dim, numComp>(meshInfo); } + /** + * @brief create a Lagrange coordinate field of the given order + * + * @details + * Each node is placed at map(x), where x is the node's position on the + * linear (straight-sided) mesh geometry. The linear field is set from the + * mesh coordinates and higher order fields are interpolated from the + * linear coordinate field, so with the default IdentityCoordinateMap the + * geometry is the mesh's own and is straight-sided. A nonlinear map gives + * curved geometry for orders above one. + * + * @tparam order the order of the coordinate field + * @param map a device callable taking and returning a + * Kokkos::Array of physical coordinates + */ + template + auto CreateLagrangeCoordinateField( + const CoordinateMap &map = CoordinateMap()) const { + auto fieldWithCtrlr = CreateLagrangeField(); + auto field = fieldWithCtrlr.field; + // note order 1 is called in the constructor to initialize the coordinate field + // it must be called before higher order coordinate fields are constructed + if constexpr (order == 1) { + const auto meshDim = meshInfo.dim; + const auto coords = mesh.coords(); + MeshField::parallel_for( + ExecutionSpace(), {0}, {meshInfo.numVtx}, + SetCoordFieldFunctor(field, coords, meshDim, map), + "setCoordField"); + } else { + const auto linear = CreateLagrangeElement(coordField.field); + auto target = CreateLagrangeElement(field, linear); + MeshField::setDofByInterpolation(target, map); + } + return fieldWithCtrlr; + } + auto getCoordField() { return coordField; } + /** + * @brief create a FieldElement over the mesh elements for a Lagrange field + * + * @details + * The shape and dof map are chosen from the field's order. + * + * @param coordinate field used to create the element + */ + template + auto CreateLagrangeElement(const Field &field) const { + static_assert(dim == 2 || dim == 3, + "CreateLagrangeElement supports triangles and tetrahedra"); + constexpr int order = Field::Order; + if constexpr (dim == 2) { + const auto [shp, map] = Omegah::getTriangleElement(mesh); + return MeshField::FieldElement(meshInfo.numTri, field, shp, map); + } else { + const auto [shp, map] = Omegah::getTetrahedronElement(mesh); + return MeshField::FieldElement(meshInfo.numTet, field, shp, map); + } + } + + /** + * @brief create a MappedElement for a field on the given geometry + * + * @param field the Lagrange field to evaluate + * @param geomElem the geometry element: CreateLagrangeElement of a + * coordinate field + */ + template + auto CreateLagrangeElement(const ShapeField &field, + const GeomElem &geomElem) const { + return MeshField::MappedElement(CreateLagrangeElement(field), geomElem); + } + // FIXME support 2d and 3d and fields with order>1 template void writeVtk(Field &field) const { using FieldDataType = typename decltype(field.vtxField)::BaseType; diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index cb2890d..6adc30c 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -157,6 +157,12 @@ struct FieldElement { const ElementDofHolderAccessor elm2dofIn) : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), elm2dof(elm2dofIn) {} + + /** @brief the topology of the element's reference element */ + static constexpr KOKKOS_FUNCTION auto getTopology() { + return ElementDofHolderAccessor::getTopology(); + } + /* general template for baseType which simply sets type */ template struct baseType { @@ -173,40 +179,144 @@ struct FieldElement { FieldAccessor::numComp>; static const size_t NumComponents = FieldAccessor::numComp; + /** @brief every dof of one element, indexed [node * NumComponents + comp] */ + using DofArray = + Kokkos::Array::type, + ShapeType::numNodes * NumComponents>; + + /** @brief parametric gradients of every component at one point, indexed + * [direction * NumComponents + comp] */ + using LocalGradArray = + Kokkos::Array::type, + MeshEntDim * NumComponents>; + /** * @brief - * evaluate the field in the specified element at the specified - * parametric/local/area coordinate - * - * @todo add expression in documetation for summation of shape function * - * field values + * gather element's dofs into an array * - * @details - * heavily based on SCOREC/core @ 7cd76473 apf/apfElement.cc * * @param ent the mesh entity index - * @param localCoord the parametric coordinate - * @return the result of evaluation + * @return the element's dofs, indexed [node * NumComponents + component] */ - KOKKOS_INLINE_FUNCTION ValArray - getValue(int ent, Kokkos::Array localCoord) const { + KOKKOS_INLINE_FUNCTION DofArray gatherDofs(int ent) const { assert(ent >= 0); assert(static_cast(ent) < numMeshEnts); - ValArray c; - const auto shapeValues = shapeFn.getValues(localCoord); - for (size_t ci = 0; ci < NumComponents; ++ci) - c[ci] = 0; + DofArray dofs; for (auto topo : elm2dof.getTopology()) { // element topology - for (size_t ni = 0; ni < shapeFn.numNodes; ++ni) { + for (size_t ni = 0; ni < ShapeType::numNodes; ++ni) { for (size_t ci = 0; ci < NumComponents; ++ci) { - auto map = elm2dof(ni, ci, ent, topo); - const auto fval = + const auto map = elm2dof(ni, ci, ent, topo); + dofs[ni * NumComponents + ci] = field(map.entity, map.node, map.component, map.topo); - c[ci] += fval * shapeValues[ni]; } } } - return c; + return dofs; + } + + /** + * @brief + * the coefficients of one element on its reference basis + * + * @details + * The coefficients depend only on the element, not on the evaluation point, + * so compute them once per element and pass them to valueFromCoefficients + * and localGradientFromCoefficients for each point. For a FieldElement the + * stored dofs are the coefficients; elements whose dofs must be transformed + * to the reference element (e.g., the Bell element) override this. + * + * @param ent the mesh entity index + * @return coefficients, indexed [node * NumComponents + component] + */ + KOKKOS_INLINE_FUNCTION DofArray coefficients(int ent) const { + return gatherDofs(ent); + } + + /** + * @brief + * contract an element's coefficients with the basis at a parametric + * coordinate + * + * @details + * heavily based on SCOREC/core @ 7cd76473 apf/apfElement.cc + * + * @param c coefficients from coefficients(ent) + * @param localCoord the parametric coordinate + * @return the field value + */ + template + KOKKOS_INLINE_FUNCTION ValArray valueFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + const auto shapeValues = shapeFn.getValues(localCoord); + ValArray v; + for (size_t ci = 0; ci < NumComponents; ++ci) + v[ci] = 0; + for (size_t ni = 0; ni < ShapeType::numNodes; ++ni) + for (size_t ci = 0; ci < NumComponents; ++ci) + v[ci] += c[ni * NumComponents + ci] * shapeValues[ni]; + return v; + } + + /** + * @brief + * contract an element's coefficients with the basis gradients at a + * parametric coordinate + * + * @details + * The result is the gradient with respect to the parametric coordinates, + * not the physical coordinates; mapping it to the physical element requires + * the geometry's Jacobian. + * + * @param c coefficients from coefficients(ent) + * @param localCoord the parametric coordinate + * @return parametric gradients, indexed [direction * NumComponents + component] + */ + template + KOKKOS_INLINE_FUNCTION LocalGradArray localGradientFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + return localGradientFromShapeGradients( + c, shapeFn.getLocalGradients(localCoord)); + } + + /** + * @brief + * localGradientFromCoefficients with the shape function gradients already + * evaluated, so they can be shared with a geometry that uses the same shape + * + * @param c coefficients from coefficients(ent) + * @param shapeGrads shapeFn.getLocalGradients at the parametric coordinate + * @return parametric gradients, indexed [direction * NumComponents + component] + */ + template + KOKKOS_INLINE_FUNCTION LocalGradArray localGradientFromShapeGradients( + Coefficients const &c, ShapeGradients const &shapeGrads) const { + LocalGradArray g; + for (size_t i = 0; i < MeshEntDim * NumComponents; ++i) + g[i] = 0; + for (size_t ni = 0; ni < ShapeType::numNodes; ++ni) + for (size_t d = 0; d < MeshEntDim; ++d) + for (size_t ci = 0; ci < NumComponents; ++ci) + g[d * NumComponents + ci] += + c[ni * NumComponents + ci] * shapeGrads[ni * MeshEntDim + d]; + return g; + } + + /** + * @brief + * evaluate the field in the specified element at the specified + * parametric/local/area coordinate + * + * @details + * Prefer coefficients and valueFromCoefficients when evaluating an element + * at more than one point. + * + * @param ent the mesh entity index + * @param localCoord the parametric coordinate + * @return the result of evaluation + */ + KOKKOS_INLINE_FUNCTION ValArray + getValue(int ent, Kokkos::Array localCoord) const { + return valueFromCoefficients(coefficients(ent), localCoord); } using NodeArray = @@ -227,6 +337,51 @@ struct FieldElement { return c; } + /** @brief the Jacobian at one point, row-major with J[i][j] = + * dx_j/dxi_i */ + using JacobianArray = Kokkos::Array; + + /** + * @brief + * compute the Jacobian at a parametric coordinate from the node coordinates + * + * @details + * The field must be a coordinate field. + * + * heavily based on SCOREC/core @ 7cd76473 apf/apfVectorElement.cc + * + * @param x node coordinates from getNodeValues + * @param localCoord the parametric coordinate + * @return the Jacobian, row-major with J[i][j] = dx_j/dxi_i + */ + KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromNodeValues( + NodeArray const &x, Kokkos::Array localCoord) const { + return jacobianFromShapeGradients(x, shapeFn.getLocalGradients(localCoord)); + } + + /** + * @brief + * jacobianFromNodeValues with the shape function gradients already + * evaluated, so they can be shared with a field that uses the same shape + * + * @param x node coordinates from getNodeValues + * @param grads shapeFn.getLocalGradients at the parametric coordinate + * @return the Jacobian, row-major with J[i][j] = dx_j/dxi_i + */ + template + KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromShapeGradients( + NodeArray const &x, ShapeGradients const &grads) const { + JacobianArray J; + for (size_t i = 0; i < MeshEntDim * MeshEntDim; ++i) + J[i] = 0; + for (size_t n = 0; n < ShapeType::numNodes; ++n) + for (size_t i = 0; i < MeshEntDim; ++i) + for (size_t j = 0; j < MeshEntDim; ++j) + J[i * MeshEntDim + j] += + grads[n * MeshEntDim + i] * x[n * MeshEntDim + j]; + return J; + } + /** * @brief * compute the Jacobian of an edge @@ -516,11 +671,12 @@ evaluate(FieldElement &fes, Kokkos::View localCoords, Kokkos::parallel_for( fes.numMeshEnts, KOKKOS_LAMBDA(const int ent) { Kokkos::Array lc; + const auto c = fes.coefficients(ent); // TODO use nested parallel for? for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) { for (size_t i = 0; i < localCoords.extent(1); i++) // better way? lc[i] = localCoords(pt, i); - const auto val = fes.getValue(ent, lc); + const auto val = fes.valueFromCoefficients(c, lc); for (size_t i = 0; i < numComponents; i++) res(pt, i) = val[i]; } @@ -599,5 +755,93 @@ Kokkos::View getJacobianDeterminants(FieldElement &fes, Matrices J) { return fes.getJacobianDeterminants(J); } +/** + * @brief + * For each dof holder of a field, the element + * that owns it; see getDofHolderOwningEntities + */ +struct DofHolderOwningEntities { + static constexpr size_t NumTopologies = Pyramid + 1; + Kokkos::Array offset; + Kokkos::View owner; + + /** @brief the index of the element that owns the given dof holder */ + KOKKOS_INLINE_FUNCTION LO operator()(Mesh_Topology topo, LO entity) const { + return owner(offset[topo] + entity); + } +}; + +/** + * @brief + * the owning entity of a dof that is shared, i.e., vertex or edge nodes + * is the shared entity with the lowest id + */ +template +DofHolderOwningEntities getDofHolderOwningEntities(FieldElement const &fes) { + constexpr size_t numNodes = decltype(fes.shapeFn)::numNodes; + constexpr size_t numTopologies = DofHolderOwningEntities::NumTopologies; + const auto &meshInfo = fes.field.meshInfo; + const LO numElements = fes.numMeshEnts; + DofHolderOwningEntities owners; + owners.offset[0] = 0; + for (size_t t = 0; t < numTopologies; ++t) + owners.offset[t + 1] = + owners.offset[t] + meshInfo.numEntities(static_cast(t)); + owners.owner = Kokkos::View( + Kokkos::view_alloc(Kokkos::WithoutInitializing, "dofHolderOwner"), + owners.offset[numTopologies]); + Kokkos::deep_copy(owners.owner, numElements); + const auto offset = owners.offset; + const auto owner = owners.owner; + Kokkos::parallel_for( + "getDofHolderOwningEntities", numElements, KOKKOS_LAMBDA(const int ent) { + for (auto topo : fes.elm2dof.getTopology()) { + for (size_t n = 0; n < numNodes; ++n) { + const auto map = fes.elm2dof(n, 0, ent, topo); + Kokkos::atomic_min(&owner(offset[map.topo] + map.entity), ent); + } + } + }); + return owners; +} + +/** + * @brief + * evaluates a source function at the target field element dof holder locations + * after this function runs, the target's field values will modified + * + * @details + * Reference space interpolation. To interpolate in physical space see MappedElement overload. + */ +template +void setDofByInterpolation(FieldElement &target, Source const &source) { + constexpr size_t numNodes = decltype(target.shapeFn)::numNodes; + constexpr size_t meshEntDim = FieldElement::MeshEntDim; + constexpr size_t numComp = FieldElement::NumComponents; + // we get a list of owners so that we only perform evaluation once + // per shared dof + const auto owners = getDofHolderOwningEntities(target); + Kokkos::parallel_for( + "interpolate", target.numMeshEnts, KOKKOS_LAMBDA(const int ent) { + const auto nodeCoords = target.shapeFn.getNodeParametricCoords(); + for (auto topo : target.elm2dof.getTopology()) { + for (size_t n = 0; n < numNodes; ++n) { + const auto holder = target.elm2dof(n, 0, ent, topo); + if (owners(holder.topo, holder.entity) != ent) + continue; + Kokkos::Array xi; + for (size_t d = 0; d < meshEntDim; ++d) + xi[d] = nodeCoords[n * meshEntDim + d]; + const auto val = source(ent, xi); + for (size_t c = 0; c < numComp; ++c) { + const auto map = target.elm2dof(n, c, ent, topo); + target.field(map.entity, map.node, map.component, map.topo) = + val[c]; + } + } + } + }); +} + } // namespace MeshField #endif diff --git a/src/MeshField_For.hpp b/src/MeshField_For.hpp index 87a981c..3b605ad 100644 --- a/src/MeshField_For.hpp +++ b/src/MeshField_For.hpp @@ -7,7 +7,7 @@ namespace MeshField { template void parallel_for(ExecutionSpace, const std::initializer_list &start, const std::initializer_list &end, - FunctorType &vectorKernel, std::string tag) { + const FunctorType &vectorKernel, std::string tag) { constexpr auto funcRank = MeshFieldUtil::function_traits::arity; assert(funcRank >= 1); Kokkos::Array a_start = diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index da80ba4..04bd759 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -253,7 +253,7 @@ class Integrator { * \param fes FieldElement * */ template void process(FieldElement &fes) { - constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology(); + constexpr auto topo = FieldElement::getTopology(); pre(); auto ip = getIntegrationPoints(order); auto localCoords = getIntegrationPointLocalCoords(fes, ip); diff --git a/src/MeshField_LinearAlgebra.hpp b/src/MeshField_LinearAlgebra.hpp new file mode 100644 index 0000000..790c10d --- /dev/null +++ b/src/MeshField_LinearAlgebra.hpp @@ -0,0 +1,72 @@ +#ifndef MESHFIELD_LINEARALGEBRA_HPP +#define MESHFIELD_LINEARALGEBRA_HPP + +#include +#include +#include + +/** + * @file MeshField_LinearAlgebra.hpp + * @brief small dense matrix operations callable on host and device + * + * Matrices are Dim x Dim and stored row-major in a + * Kokkos::Array. + */ + +namespace MeshField { + +/** + * @brief the determinant of a 1x1, 2x2, or 3x3 matrix stored row-major + */ +template +KOKKOS_INLINE_FUNCTION Real +determinant(Kokkos::Array const &A) { + static_assert(Dim >= 1 && Dim <= 3, "only 1x1, 2x2, and 3x3 are supported"); + if constexpr (Dim == 1) { + return A[0]; + } else if constexpr (Dim == 2) { + return A[0] * A[3] - A[1] * A[2]; + } else { + return A[0] * (A[4] * A[8] - A[5] * A[7]) - + A[1] * (A[3] * A[8] - A[5] * A[6]) + + A[2] * (A[3] * A[7] - A[4] * A[6]); + } +} + +/** + * @brief invert a 1x1, 2x2, or 3x3 matrix stored row-major + */ +template +KOKKOS_INLINE_FUNCTION Kokkos::Array +invert(Kokkos::Array const &A) { + static_assert(Dim >= 1 && Dim <= 3, "only 1x1, 2x2, and 3x3 are supported"); + const auto det = determinant(A); + assert(Kokkos::fabs(det) > Epsilon); + Kokkos::Array inv; + if constexpr (Dim == 1) { + inv[0] = 1; + } else if constexpr (Dim == 2) { + inv[0] = A[3]; + inv[1] = -A[1]; + inv[2] = -A[2]; + inv[3] = A[0]; + } else { + // adjugate + inv[0] = A[4] * A[8] - A[5] * A[7]; + inv[1] = A[2] * A[7] - A[1] * A[8]; + inv[2] = A[1] * A[5] - A[2] * A[4]; + inv[3] = A[5] * A[6] - A[3] * A[8]; + inv[4] = A[0] * A[8] - A[2] * A[6]; + inv[5] = A[2] * A[3] - A[0] * A[5]; + inv[6] = A[3] * A[7] - A[4] * A[6]; + inv[7] = A[1] * A[6] - A[0] * A[7]; + inv[8] = A[0] * A[4] - A[1] * A[3]; + } + for (size_t i = 0; i < Dim * Dim; ++i) + inv[i] /= det; + return inv; +} + +} // namespace MeshField + +#endif diff --git a/src/MeshField_MappedElement.hpp b/src/MeshField_MappedElement.hpp new file mode 100644 index 0000000..93d1d9f --- /dev/null +++ b/src/MeshField_MappedElement.hpp @@ -0,0 +1,269 @@ +#ifndef MESHFIELD_MAPPEDELEMENT_HPP +#define MESHFIELD_MAPPEDELEMENT_HPP + +#include +#include +#include +#include +#include +#include +#include // std::declval + +/** + * @file MeshField_MappedElement.hpp + * @brief Elements whose geometry, dofs, or values need mapping between the + * reference element and the physical element. + * + * A MappedElement composes three things: + * + * - **a geometry element**: a FieldElement over the coordinate field that + * supplies the map from the reference element to the physical element. + * - **a DofTransformation**: converts the dofs stored on mesh entities to + * coefficients on the reference basis. See: + * R.C. Kirby, "A general approach to transforming finite + * elements", SMAI J. Comput. Math. 4 (2018) 197-224 , and https://defelement.org/ + * - **a MapType**: converts values and derivatives on the reference element + * to the physical element. See: https://defelement.org + */ + +namespace MeshField { + +struct IdentityDofTransformation { + template + KOKKOS_INLINE_FUNCTION auto coefficients(FieldElem const &field, + GeomElem const &, int ent) const { + return field.coefficients(ent); + } +}; + +struct IdentityMap { + template + KOKKOS_INLINE_FUNCTION ValArray value(ValArray const &ref) const { + return ref; + } + + /** + * @brief computes the physical gradient + * + * @details + * With J[i][j] = dx_j/dxi_i the chain rule gives grad_xi = J grad_x, so + * grad_x = J^{-1} grad_xi + */ + template + KOKKOS_INLINE_FUNCTION GradArray gradient(GradArray const &refGrad, + JacobianArray const &J) const { + // sqrt not constexpr until C++26. Is there a more sane way to do this? + constexpr size_t Dim = JacobianArray::size() == 1 ? 1 + : JacobianArray::size() == 4 ? 2 + : 3; + static_assert(Dim * Dim == JacobianArray::size()); + constexpr auto numComp = GradArray::size() / Dim; + const auto Jinv = invert(J); + GradArray g; + for (size_t i = 0; i < Dim; ++i) { + for (size_t ci = 0; ci < numComp; ++ci) { + Real sum = 0; + for (size_t j = 0; j < Dim; ++j) + sum += Jinv[i * Dim + j] * refGrad[j * numComp + ci]; + g[i * numComp + ci] = sum; + } + } + return g; + } +}; + +/** + * @brief Element type used for non-isoparametric, and other elements that require mapping + * + * @details + * The field element supplies the basis and the stored dofs. The geometry + * element supplies the map from the reference element to the physical + * element. I + * + * Evaluation: + * @code + * const auto c = elm.coefficients(ent); // once per element + * const auto g = elm.geometry(ent); // if needed + * for (pt : points) { + * elm.valueFromCoefficients(c, xi); // values + * elm.gradientFromCoefficients(c, g, xi); // gradients + * } + * @endcode + * + * @tparam FieldElem a FieldElement over the field being evaluated + * @tparam GeomElem a FieldElement over the coordinate field + * @tparam DofTransformation see IdentityDofTransformation + * @tparam MapType see IdentityMap + */ +template +struct MappedElement { + static_assert(FieldElem::MeshEntDim == GeomElem::MeshEntDim, + "the field and geometry elements must have the same " + "parametric dimension"); + static_assert(FieldElem::getTopology()[0] == GeomElem::getTopology()[0], + "the field and geometry elements must have the same topology"); + + FieldElem field; + GeomElem geom; + DofTransformation dofTransformation; + MapType mapType; + const size_t numMeshEnts; + + static constexpr size_t MeshEntDim = FieldElem::MeshEntDim; + static constexpr size_t NumComponents = FieldElem::NumComponents; + using ValArray = typename FieldElem::ValArray; + + /** @brief the topology of the reference element */ + static constexpr KOKKOS_FUNCTION auto getTopology() { + return GeomElem::getTopology(); + } + + /** @brief physical gradients + * [direction * NumComponents + component] */ + using GradArray = typename FieldElem::LocalGradArray; + /** @brief the Jacobian, J[i][j] = dx_j/dxi_i */ + using JacobianArray = typename GeomElem::JacobianArray; + + /** @brief the geometry map is affine, so the Jacobian is constant over + * each element */ + static constexpr bool IsAffine = decltype(GeomElem::shapeFn)::Order == 1; + /** @brief the field and geometry use the same shape functions and dof map */ + static constexpr bool IsIsoparametric = + std::is_same_v && + std::is_same_v; + + /** @brief the field's coefficients on the reference basis */ + using Coefficients = decltype(std::declval().coefficients( + std::declval(), std::declval(), 0)); + using Geometry = + std::conditional_t; + + MappedElement( + FieldElem const &fieldIn, GeomElem const &geomIn, + DofTransformation const &dofTransformationIn = DofTransformation(), + MapType const &mapTypeIn = MapType()) + : field(fieldIn), geom(geomIn), dofTransformation(dofTransformationIn), + mapType(mapTypeIn), numMeshEnts(fieldIn.numMeshEnts) { + if (geomIn.numMeshEnts != fieldIn.numMeshEnts) { + fail("the field and geometry elements have the same number of " + "mesh entities (%zu != %zu)\n", + fieldIn.numMeshEnts, geomIn.numMeshEnts); + } + } + + KOKKOS_INLINE_FUNCTION Coefficients coefficients(int ent) const { + return dofTransformation.coefficients(field, geom, ent); + } + + /** + * @brief the geometry of one element + */ + KOKKOS_INLINE_FUNCTION Geometry geometry(int ent) const { + if constexpr (IsAffine) { + return geom.jacobianFromNodeValues(geom.getNodeValues(ent), + Kokkos::Array{}); + } else { + return geom.getNodeValues(ent); + } + } + + /** @brief the physical value of the field at a parametric coordinate */ + KOKKOS_INLINE_FUNCTION ValArray valueFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + return mapType.value(field.valueFromCoefficients(c, localCoord)); + } + + /** @brief the Jacobian of the geometry at a parametric coordinate */ + KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromGeometry( + Geometry const &g, Kokkos::Array localCoord) const { + if constexpr (IsAffine) { + return g; + } else { + return geom.jacobianFromNodeValues(g, localCoord); + } + } + + /** @brief the physical gradient of the field at a parametric coordinate */ + KOKKOS_INLINE_FUNCTION GradArray + gradientFromCoefficients(Coefficients const &c, Geometry const &g, + Kokkos::Array localCoord) const { + if constexpr (IsIsoparametric && !IsAffine) { + // the field and the geometry share the shape function gradients + const auto grads = field.shapeFn.getLocalGradients(localCoord); + return mapType.gradient(field.localGradientFromShapeGradients(c, grads), + geom.jacobianFromShapeGradients(g, grads)); + } else { + return mapType.gradient( + field.localGradientFromCoefficients(c, localCoord), + jacobianFromGeometry(g, localCoord)); + } + } + + /** + * @brief evaluate the field at a single parametric coordinate + * + * Prefer coefficients and valueFromCoefficients when evaluating an element + * at more than one point. + */ + KOKKOS_INLINE_FUNCTION ValArray + getValue(int ent, Kokkos::Array localCoord) const { + return valueFromCoefficients(coefficients(ent), localCoord); + } + + /** @brief Jacobians, from the geometry element */ + Kokkos::View getJacobians(Kokkos::View localCoords, + Kokkos::View offsets) const { + return geom.getJacobians(localCoords, offsets); + } + + /** @brief Jacobian determinants, from the geometry element */ + template + Kokkos::View getJacobianDeterminants(Matrices const &J) const { + return geom.getJacobianDeterminants(J); + } +}; + +/** + * @brief + * evaluates a source function at the dof holder locations of the target's + * field, in reference space + */ +template + requires std::is_invocable_v const &> +void setDofByInterpolation( + MappedElement &target, + Source const &source) { + setDofByInterpolation(target.field, source); +} + +/** + * @brief + * evaluates a function at the dof holder locations of the target's field, in + * physical space + * + */ +template + requires std::is_invocable_v const &> +void setDofByInterpolation( + MappedElement &target, + Function const &f) { + const auto geom = target.geom; + setDofByInterpolation( + target.field, + KOKKOS_LAMBDA(const int ent, + Kokkos::Array const &xi) { + return f(geom.getValue(ent, xi)); + }); +} + +} // namespace MeshField + +#endif diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 5175b02..6178a80 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -17,15 +17,38 @@ namespace MeshField { * On-process mesh metadata */ struct MeshInfo { - int numVtx; // entDim = 0 - int numEdge; // entDim = 1 - int numTri; // entDim = 2 - int numQuad; // entDim = 2 - int numTet; // entDim = 3 - int numHex; // entDim = 3 - int numPrism; // entDim = 3 - int numPyramid; // entDim = 3 - int dim; + int numVtx = 0; // entDim = 0 + int numEdge = 0; // entDim = 1 + int numTri = 0; // entDim = 2 + int numQuad = 0; // entDim = 2 + int numTet = 0; // entDim = 3 + int numHex = 0; // entDim = 3 + int numPrism = 0; // entDim = 3 + int numPyramid = 0; // entDim = 3 + int dim = 0; + + /** @brief runtime map of number of entities based on topology */ + int numEntities(Mesh_Topology topo) const { + switch (topo) { + case Vertex: + return numVtx; + case Edge: + return numEdge; + case Triangle: + return numTri; + case Quad: + return numQuad; + case Tetrahedron: + return numTet; + case Hexahedron: + return numHex; + case Prism: + return numPrism; + case Pyramid: + return numPyramid; + } + return 0; + } }; /** diff --git a/test/MeshField_TestFields.hpp b/test/MeshField_TestFields.hpp new file mode 100644 index 0000000..4f54fa6 --- /dev/null +++ b/test/MeshField_TestFields.hpp @@ -0,0 +1,155 @@ +#ifndef MESHFIELD_TESTFIELDS_HPP +#define MESHFIELD_TESTFIELDS_HPP + +/** + * @file MeshField_TestFields.hpp + * @brief controllers with one field of each rank, and helpers that fill and + * check a Field or a Kokkos::View of any rank, shared by the controller tests + */ + +#include "MeshField_TestHarness.hpp" +#include +#include +#include + +using MemorySpace = ExecutionSpace::memory_space; + +/** @brief the extents of a field or View of rank up to 5 */ +using Extents = Kokkos::Array; + +// Fields of type T, one of each rank up to the controller's MAX_RANK: 5 for +// Kokkos, 4 for Cabana. Each field has a different number of entries (its +// first extent), so mixed up fields show. The Kokkos data types include the +// first extent, Cabana's do not. +constexpr int Y = 3, Z = 4, A = 5, B = 2; +constexpr Extents rank1 = {10}, rank2 = {11, Y}, rank3 = {12, Y, Z}, + rank4 = {13, Y, Z, A}, rank5 = {14, Y, Z, A, B}; +template +using KokkosRanks = + MeshField::KokkosController; +#ifdef MESHFIELDS_ENABLE_CABANA +template +using CabanaRanks = MeshField::CabanaController; +#endif + +/** @brief create a controller with the fields of ranks 1 to MAX_RANK */ +template Ctrlr createRanks() { + // both controllers take the runtime extents, here only the first ones + if constexpr (Ctrlr::MAX_RANK >= 5) + return Ctrlr({rank1[0], rank2[0], rank3[0], rank4[0], rank5[0]}); + else + return Ctrlr({rank1[0], rank2[0], rank3[0], rank4[0]}); +} + +/** @brief a value that differs at every index, so transposed indices show; + * every value used is exact in float and double */ +KOKKOS_INLINE_FUNCTION int indexValue(int i, int j = 0, int k = 0, int l = 0, + int m = 0) { + return i + 10 * j + 100 * k + 1000 * l + 10000 * m; +} + +/** @brief set every entry of a Field or View of rank R to indexValue + + * offset */ +template +void fill(Container const &c, Extents n, int offset = 0) { + static_assert(R >= 1 && R <= 5); + if constexpr (R == 1) { + auto kernel = KOKKOS_LAMBDA(int i) { c(i) = indexValue(i) + offset; }; + MeshField::parallel_for(ExecutionSpace(), {0}, {n[0]}, kernel, "fill1"); + } else if constexpr (R == 2) { + auto kernel = KOKKOS_LAMBDA(int i, int j) { + c(i, j) = indexValue(i, j) + offset; + }; + MeshField::parallel_for(ExecutionSpace(), {0, 0}, {n[0], n[1]}, kernel, + "fill2"); + } else if constexpr (R == 3) { + auto kernel = KOKKOS_LAMBDA(int i, int j, int k) { + c(i, j, k) = indexValue(i, j, k) + offset; + }; + MeshField::parallel_for(ExecutionSpace(), {0, 0, 0}, {n[0], n[1], n[2]}, + kernel, "fill3"); + } else if constexpr (R == 4) { + auto kernel = KOKKOS_LAMBDA(int i, int j, int k, int l) { + c(i, j, k, l) = indexValue(i, j, k, l) + offset; + }; + MeshField::parallel_for(ExecutionSpace(), {0, 0, 0, 0}, + {n[0], n[1], n[2], n[3]}, kernel, "fill4"); + } else { + auto kernel = KOKKOS_LAMBDA(int i, int j, int k, int l, int m) { + c(i, j, k, l, m) = indexValue(i, j, k, l, m) + offset; + }; + MeshField::parallel_for(ExecutionSpace(), {0, 0, 0, 0, 0}, + {n[0], n[1], n[2], n[3], n[4]}, kernel, "fill5"); + } +} + +/** @brief the number of entries of a Field or View of rank R that differ from + * indexValue */ +template +MeshField::LO mismatches(Container const &c, Extents n) { + static_assert(R >= 1 && R <= 5); + MeshField::LO errors = 0; + if constexpr (R == 1) { + auto kernel = KOKKOS_LAMBDA(int i, MeshField::LO &lerr) { + lerr += c(i) != indexValue(i); + }; + MeshField::parallel_reduce(ExecutionSpace(), "mismatches1", {0}, {n[0]}, + kernel, errors); + } else if constexpr (R == 2) { + auto kernel = KOKKOS_LAMBDA(int i, int j, MeshField::LO &lerr) { + lerr += c(i, j) != indexValue(i, j); + }; + MeshField::parallel_reduce(ExecutionSpace(), "mismatches2", {0, 0}, + {n[0], n[1]}, kernel, errors); + } else if constexpr (R == 3) { + auto kernel = KOKKOS_LAMBDA(int i, int j, int k, MeshField::LO &lerr) { + lerr += c(i, j, k) != indexValue(i, j, k); + }; + MeshField::parallel_reduce(ExecutionSpace(), "mismatches3", {0, 0, 0}, + {n[0], n[1], n[2]}, kernel, errors); + } else if constexpr (R == 4) { + auto kernel = + KOKKOS_LAMBDA(int i, int j, int k, int l, MeshField::LO &lerr) { + lerr += c(i, j, k, l) != indexValue(i, j, k, l); + }; + MeshField::parallel_reduce(ExecutionSpace(), "mismatches4", {0, 0, 0, 0}, + {n[0], n[1], n[2], n[3]}, kernel, errors); + } else { + auto kernel = KOKKOS_LAMBDA(int i, int j, int k, int l, int m, + MeshField::LO &lerr) { + lerr += c(i, j, k, l, m) != indexValue(i, j, k, l, m); + }; + MeshField::parallel_reduce(ExecutionSpace(), "mismatches5", + {0, 0, 0, 0, 0}, {n[0], n[1], n[2], n[3], n[4]}, + kernel, errors); + } + return errors; +} + +/** @brief a Kokkos::View of rank R */ +template auto createView(Extents n) { + static_assert(R >= 1 && R <= 5); + if constexpr (R == 1) + return Kokkos::View("view1", n[0]); + else if constexpr (R == 2) + return Kokkos::View("view2", n[0], n[1]); + else if constexpr (R == 3) + return Kokkos::View("view3", n[0], n[1], n[2]); + else if constexpr (R == 4) + return Kokkos::View("view4", n[0], n[1], n[2], n[3]); + else + return Kokkos::View("view5", n[0], n[1], n[2], n[3], n[4]); +} + +/** @brief the number of extents of a field of rank R that differ from n */ +template +MeshField::LO extentMismatches(Field const &field, Extents n) { + MeshField::LO errors = 0; + for (int d = 0; d < R; ++d) + errors += static_cast(field.size(d)) != n[d]; + return errors; +} + +#endif diff --git a/test/MeshField_TestFunctions.hpp b/test/MeshField_TestFunctions.hpp new file mode 100644 index 0000000..4b75d4f --- /dev/null +++ b/test/MeshField_TestFunctions.hpp @@ -0,0 +1,129 @@ +#ifndef MESHFIELD_TESTFUNCTIONS_HPP +#define MESHFIELD_TESTFUNCTIONS_HPP + +/** + * @file MeshField_TestFunctions.hpp + * @brief analytic functions and coordinate maps shared by the tests + */ + +#include +#include + +template using Vec = Kokkos::Array; + +/** @brief f(x) = 1 + sum_d (d+1) x_d */ +template struct LinearFunction { + static constexpr size_t degree = 1; + KOKKOS_INLINE_FUNCTION MeshField::Real operator()(Vec const &x) const { + MeshField::Real f = 1; + for (size_t d = 0; d < dim; ++d) + f += (d + 1) * x[d]; + return f; + } + KOKKOS_INLINE_FUNCTION Vec grad(Vec const &) const { + Vec g; + for (size_t d = 0; d < dim; ++d) + g[d] = d + 1; + return g; + } +}; + +/** + * @brief f(x) = 1 + sum_d ((d+1) x_d + x_d^2) + x_0 x_1 + * + * The cross term gives the Hessian an off-diagonal entry. + */ +template struct QuadraticFunction { + static_assert(dim >= 2, "the cross term needs at least two dimensions"); + static constexpr size_t degree = 2; + KOKKOS_INLINE_FUNCTION MeshField::Real operator()(Vec const &x) const { + MeshField::Real f = 1 + x[0] * x[1]; + for (size_t d = 0; d < dim; ++d) + f += (d + 1) * x[d] + x[d] * x[d]; + return f; + } + KOKKOS_INLINE_FUNCTION Vec grad(Vec const &x) const { + Vec g; + for (size_t d = 0; d < dim; ++d) + g[d] = (d + 1) + 2 * x[d]; + g[0] += x[1]; + g[1] += x[0]; + return g; + } +}; + +/** + * @brief the value of a numComp component function whose component c is + * (c+1)*f, which checks the component ordering + */ +template +KOKKOS_INLINE_FUNCTION Kokkos::Array +valueScaledByComponent(Function const &f, Vec const &x) { + Kokkos::Array val; + for (size_t c = 0; c < numComp; ++c) + val[c] = (c + 1) * f(x); + return val; +} + +/** + * @brief the gradient of the function in valueScaledByComponent, indexed + * [direction * numComp + component]; f must provide grad + */ +template +KOKKOS_INLINE_FUNCTION Kokkos::Array +gradScaledByComponent(Function const &f, Vec const &x) { + const auto gf = f.grad(x); + Kokkos::Array g; + for (size_t d = 0; d < dim; ++d) + for (size_t c = 0; c < numComp; ++c) + g[d * numComp + c] = (c + 1) * gf[d]; + return g; +} + +/** + * @brief CoordinateMap for an affine shear that makes every Jacobian a full + * matrix: x += 0.3y + 0.1z, y += 0.2z + * + * Its determinant is one, so it preserves measure. Coordinate fields of any + * order reproduce it exactly. + */ +template struct Shear { + KOKKOS_INLINE_FUNCTION Vec operator()(Vec x) const { + x[0] += 0.3 * x[1]; + if constexpr (dim == 3) { + x[0] += 0.1 * x[2]; + x[1] += 0.2 * x[2]; + } + return x; + } +}; + +/** + * @brief CoordinateMap that bends the mesh: y += 0.2x^2 + * + * Quadratic, so quadratic coordinate fields reproduce it exactly. + */ +template struct Bend { + KOKKOS_INLINE_FUNCTION Vec operator()(Vec x) const { + x[1] += 0.2 * x[0] * x[0]; + return x; + } +}; + +/** + * @brief CoordinateMap that bulges the top of the unit square into + * y = 1 + delta x(1 - x): y *= 1 + delta x(1 - x) + * + * The top edges are parabolas in their edge parameter, so quadratic + * coordinate fields reproduce the boundary, and the area, 1 + delta/6, + * exactly. + */ +template struct Bulge { + MeshField::Real delta = 0.1; + KOKKOS_INLINE_FUNCTION Vec operator()(Vec x) const { + x[1] *= 1 + delta * x[0] * (1 - x[0]); + return x; + } +}; + +#endif diff --git a/test/MeshField_TestHarness.hpp b/test/MeshField_TestHarness.hpp new file mode 100644 index 0000000..5a99934 --- /dev/null +++ b/test/MeshField_TestHarness.hpp @@ -0,0 +1,162 @@ +#ifndef MESHFIELD_TESTHARNESS_HPP +#define MESHFIELD_TESTHARNESS_HPP + +/** + * @file MeshField_TestHarness.hpp + * @brief a test runner, reporting, and integrators shared by the tests + * + * A test's main creates a TestRunner and loops over controllers, meshes, and + * types (e.g., element factories) with its forEach functions. The loop + * bodies are template lambdas: + * + * @code + * TestRunner t(argc, argv); + * t.forEachController([&]