diff --git a/CMakeLists.txt b/CMakeLists.txt index f295708..b385901 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,7 @@ project(meshfields VERSION 1.0.1 LANGUAGES CXX) include(GNUInstallDirs) include(CMakePackageConfigHelpers) +include(CheckGit.cmake) option(MeshFields_USE_Cabana "Build with the Cabana storage backend" OFF) @@ -76,6 +77,9 @@ configure_file( "${CMAKE_CURRENT_SOURCE_DIR}/src/MeshField_Config.hpp.in" "${CMAKE_CURRENT_BINARY_DIR}/MeshField_Config.hpp") + +CheckGitSetup() +target_sources(meshfields PRIVATE ${CMAKE_CURRENT_BINARY_DIR}/MeshField_GetHash.cpp) if(MeshFields_USE_Cabana) target_link_libraries(meshfields INTERFACE Cabana::Core) target_compile_definitions(meshfields INTERFACE MESHFIELDS_ENABLE_CABANA) @@ -212,6 +216,9 @@ meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp) meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp) meshfields_add_exe(ExceptionTest test/testExceptions.cpp) meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) +meshfields_add_exe(IntegratorTest test/testIntegrator.cpp) +meshfields_add_exe(BoxIntegratorTest test/testBoxIntegrator.cpp) + meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) meshfields_add_exe(Interpolate test/testInterpolate.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) @@ -220,6 +227,7 @@ 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) + meshfields_add_exe(IntegratorPerformance test/testIntegratorPerformance.cpp) smoke_test_func(CabanaTests ./CabanaTests) test_func(ControllerPerformance ./ControllerPerformance) test_func(CabanaControllerLifetimeTest ./CabanaControllerLifetimeTest) @@ -233,6 +241,8 @@ test_func(ElementJacobian1d ./ElementJacobian1d) test_func(ElementJacobian2d ./ElementJacobian2d) test_func(ElementJacobian3d ./ElementJacobian3d) test_func(CountIntegrator ./CountIntegrator) +test_func(IntegratorTest ./IntegratorTest) +test_func(BoxIntegratorTest ./BoxIntegratorTest) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) @@ -244,7 +254,7 @@ else() will_fail_test_func(ExceptionTest ./ExceptionTest) endif() -meshfields_add_exe(OmegahSprAdapt test/testSprThwaitesAdapt.cpp) +meshfields_add_exe(OmegahSprAdapt test/testSprThwaitesAdapt.cpp) test_func(OmegahSprAdapt ./OmegahSprAdapt ${CMAKE_SOURCE_DIR}/meshes/thwaites_basal_effectiveStrain.osh/ thwaitesAdapted #output prefix diff --git a/CheckGit.cmake b/CheckGit.cmake new file mode 100644 index 0000000..a558c57 --- /dev/null +++ b/CheckGit.cmake @@ -0,0 +1,61 @@ +set(pre_configure_file ${CMAKE_CURRENT_SOURCE_DIR}/src/MeshField_GetHash.cpp.in) +set(post_configure_file ${CMAKE_CURRENT_BINARY_DIR}/MeshField_GetHash.cpp) + +function(CheckGitWrite git_hash) + file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/git-state.txt "${git_hash}\n") +endfunction() + +function(CheckGitRead git_hash) + if (EXISTS ${CMAKE_CURRENT_BINARY_DIR}/git-state.txt) + file(STRINGS ${CMAKE_CURRENT_BINARY_DIR}/git-state.txt CONTENT) + LIST(GET CONTENT 0 var) + + set(${git_hash} ${var} PARENT_SCOPE) + endif () +endfunction() + +function(CheckGitVersion) + execute_process(COMMAND git rev-parse HEAD + WORKING_DIRECTORY ${CMAKE_CURRENT_LIST_DIR} + OUTPUT_VARIABLE GIT_HASH + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + CheckGitRead(GIT_HASH_CACHE) + if (NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}) + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + endif () + + if (NOT EXISTS ${CMAKE_CURRENT_BINARY_DIR}/MeshField_GetHash.hpp) + file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/src/MeshField_GetHash.hpp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + endif() + + if (NOT DEFINED GIT_HASH_CACHE) + set(GIT_HASH_CACHE "INVALID") + endif () + + if (NOT ${GIT_HASH} STREQUAL ${GIT_HASH_CACHE} OR NOT EXISTS ${post_configure_file}) + CheckGitWrite(${GIT_HASH}) + + configure_file(${pre_configure_file} ${post_configure_file} @ONLY) + endif () + +endfunction() + +function(CheckGitSetup) + + add_custom_target(AlwaysCheckGit COMMAND ${CMAKE_COMMAND} + -DRUN_CHECK_GIT_VERSION=1 + -Dpre_configure_dir=${CMAKE_CURRENT_SOURCE_DIR} + -Dpost_configure_file=${CMAKE_CURRENT_BINARY_DIR} + -DGIT_HASH_CACHE=${GIT_HASH_CACHE} + -P ${CMAKE_CURRENT_SOURCE_DIR}/CheckGit.cmake + BYPRODUCTS ${post_configure_file} + ) + + + CheckGitVersion() +endfunction() + +if (RUN_CHECK_GIT_VERSION) + CheckGitVersion() +endif () diff --git a/docs/Ordering.tex b/docs/Ordering.tex new file mode 100644 index 0000000..89e536d --- /dev/null +++ b/docs/Ordering.tex @@ -0,0 +1,87 @@ +\documentclass{article} +\usepackage[letterpaper,textwidth=5in,right=2.5in,textheight=9in]{geometry} +\usepackage{tikz} +\pagestyle{empty} +\begin{document} +\noindent +Joshua Kloepfer\\ +9/30/2025\\ +\textbf{Meshfields DOF Holder and Integration Point Ordering}\\\\ +The node ordering and parametric coordinates below follow +Zienkiewicz, Taylor, and Zhu, ``The Finite Element Method: Its Basis and +Fundamentals'', 2013, and match the shape functions defined in +\texttt{src/MeshField\_Shape.hpp}.\\ +\\ +Linear Triangle:\\ +\begin{tikzpicture}[scale=0.8] + \draw[gray, thick] (0, 0) -- (0, 5); + \draw[gray, thick] (0, 0) -- (5, 0); + \draw[gray, thick] (0, 5) -- (5, 0); + \draw[->] (-1, 0) -- (-1, 2) node[anchor=east]{$\xi_1$}; + \draw[->] (0, -1) -- (2, -1) node[anchor=north]{$\xi_0$}; + \filldraw[black] (0, 0) circle (2pt) node[anchor=north east]{V0 $(0,0)$}; + \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1 $(1,0)$}; + \filldraw[black] (0, 5) circle (2pt) node[anchor=south]{V2 $(0,1)$}; +\end{tikzpicture} +\\ +Quadratic Triangle\\ +\begin{tikzpicture}[scale=0.8] + \draw[gray, thick] (0, 0) -- (0, 5); + \draw[gray, thick] (0, 0) -- (5, 0); + \draw[gray, thick] (0, 5) -- (5, 0); + \draw[->] (-1, 0) -- (-1, 2) node[anchor=east]{$\xi_1$}; + \draw[->] (0, -1) -- (2, -1) node[anchor=north]{$\xi_0$}; + \filldraw[black] (0, 0) circle (2pt) node[anchor=north east]{V0 $(0,0)$}; + \filldraw[black] (5, 0) circle (2pt) node[anchor=west]{V1 $(1,0)$}; + \filldraw[black] (0, 5) circle (2pt) node[anchor=south]{V2 $(0,1)$}; + \filldraw[black] (2.5, 0) circle (2pt) node[anchor=north]{E0}; + \filldraw[black] (2.5, 2.5) circle (2pt) node[anchor=south west]{E1}; + \filldraw[black] (0, 2.5) circle (2pt) node[anchor=east]{E2}; +\end{tikzpicture} +\\ +Edge nodes: E0 on V0--V1, E1 on V1--V2, E2 on V0--V2; +these are nodes 3, 4 and 5 of \texttt{QuadraticTriangleShape}.\\ +Linear Tetrahedron:\\ +\begin{tikzpicture}[scale=0.8] + \draw[gray, thick] (0, 0, 0) -- (5, 0, 0); + \draw[gray, thick] (0, 0, 0) -- (0, 5, 0); + \draw[gray, thick] (5, 0, 0) -- (0, 5, 0); + \draw[gray, thick] (0, 0, 0) -- (0, 0, 5); + \draw[gray, thick] (5, 0, 0) -- (0, 0, 5); + \draw[gray, thick] (0, 5, 0) -- (0, 0, 5); + \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=north east]{V0 $(0,0,0)$}; + \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1 $(1,0,0)$}; + \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V2 $(0,1,0)$}; + \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V3 $(0,0,1)$}; +\end{tikzpicture} +\\ +Vertex node $i$ is at $L_i = 1$: V0 at the origin, V1 on $\xi_0$, +V2 on $\xi_1$, V3 on $\xi_2$. Vertex ordering must follow the right hand +rule, so that the Jacobian determinant is positive.\\ +Quadratic Tetrahedron:\\ +\begin{tikzpicture}[scale=0.8] + \draw[gray, thick] (0, 0, 0) -- (5, 0, 0); + \draw[gray, thick] (0, 0, 0) -- (0, 5, 0); + \draw[gray, thick] (5, 0, 0) -- (0, 5, 0); + \draw[gray, thick] (0, 0, 0) -- (0, 0, 5); + \draw[gray, thick] (5, 0, 0) -- (0, 0, 5); + \draw[gray, thick] (0, 5, 0) -- (0, 0, 5); + \filldraw[black] (0, 0, 0) circle (2pt) node[anchor=north east]{V0}; + \filldraw[black] (5, 0, 0) circle (2pt) node[anchor=west]{V1}; + \filldraw[black] (0, 5, 0) circle (2pt) node[anchor=south]{V2}; + \filldraw[black] (0, 0, 5) circle (2pt) node[anchor=east]{V3}; + + \filldraw[black] (2.5, 0, 0) circle (2pt) node[anchor=north]{E0}; + \filldraw[black] (0, 2.5, 0) circle (2pt) node[anchor=west]{E1}; + \filldraw[black] (0, 0, 2.5) circle (2pt) node[anchor=east]{E2}; + \filldraw[black] (2.5, 2.5, 0) circle (2pt) node[anchor=west]{E3}; + \filldraw[black] (0, 2.5, 2.5) circle (2pt) node[anchor=west]{E4}; + \filldraw[black] (2.5, 0, 2.5) circle (2pt) node[anchor=north]{E5}; +\end{tikzpicture} +\\ +Edge nodes, corresponding to nodes 4--9 of +\texttt{QuadraticTetrahedronShape}: E0 on V0--V1, E1 on V0--V2, +E2 on V0--V3, E3 on V1--V2, E4 on V2--V3, E5 on V1--V3. +Vertex ordering must follow the right hand rule and edge ordering must +follow the diagram. +\end{document} diff --git a/docs/addingShapeFunctionsForOmegah.md b/docs/addingShapeFunctionsForOmegah.md index ff3f2bc..638abea 100644 --- a/docs/addingShapeFunctionsForOmegah.md +++ b/docs/addingShapeFunctionsForOmegah.md @@ -52,6 +52,9 @@ omitted: | Triangle (2D) | \f$(\xi_0,\xi_1)\f$ | \f$[0,1]^2\f$, \f$\xi_0+\xi_1 \le 1\f$ | | Tet (3D) | \f$(\xi_0,\xi_1,\xi_2)\f$ | \f$[0,1]^3\f$, \f$\sum \xi_i \le 1\f$ | +Diagrams of the vertex and edge node ordering, with parametric coordinates, for +the linear and quadratic triangle and tetrahedron are in `docs/Ordering.tex`. + --- ## Step 1 - Define the Shape Function Struct {#step1} @@ -107,10 +110,8 @@ The mapping struct lives in `src/MeshField.hpp` inside (as numbered by the shape function) into the on-process mesh entity indices (as numbered by Omega\_h) that hold the corresponding DOFs. -### Vertex-Ordering Correction +### Mapping Omega\_h Vertices and Edges to Nodes -As previously stated, Omega\_h numbers vertices and edges within an element -differently from the MeshFields canonical ordering used by the shape functions. The Omega\_h canonical orderings for simplices are depicted below: Omega_h triangle ordering @@ -124,27 +125,50 @@ Omega_h tetrahedron ordering. Meshes composed of hypercubes (quadrilaterals and hexahedrons), pyramids, and prisms are supported by Omega_h, but they cannot be adapted. Version 1.1.0 of MeshFields only supports Omega\_h simplices. -The existing linear-triangle and linear-tetrahedron mappings correct for -the difference in Omega\_h vs MeshFields ordering with a cyclic rotation. +**Vertices.** The MeshFields vertex node ordering for triangles and tetrahedra +matches the Omega\_h element-to-vertex ordering, so vertex node \f$i\f$ maps +directly to local vertex \f$i\f$ with no rotation: ```cpp -// For triangles (triDim=2, vtxDim=0): +// For triangles (triDim=2) and tetrahedra (tetDim=3), vtxDim=0: const auto localVtxIdx = - (Omega_h::simplex_down_template(triDim, vtxDim, nodeIdx, /*ignored=*/-1) + 2) % 3; + Omega_h::simplex_down_template(elemDim, vtxDim, nodeIdx, /*ignored=*/-1); +``` -// For tetrahedra (tetDim=3, vtxDim=0): -const auto localVtxIdx = - (Omega_h::simplex_down_template(tetDim, vtxDim, nodeIdx, /*ignored=*/-1) + 3) % 4; +**Edges.** Edge nodes do not necessarily follow the Omega\_h element-to-edge +ordering, so the mapping struct uses a lookup table from node index to the +Omega\_h local edge index. +For the quadratic triangle the MeshFields edges (E0: V0-V1, E1: V1-V2, +E2: V0-V2) coincide with the Omega\_h edges, so the table is the identity. +For the quadratic tetrahedron the MeshFields edges (nodes 4-9: V0-V1, V0-V2, +V0-V3, V1-V2, V2-V3, V1-V3) are reordered to the Omega\_h edges +(0: V0-V1, 1: V1-V2, 2: V2-V0, 3: V0-V3, 4: V1-V3, 5: V2-V3): + +```cpp +// from QuadraticTetrahedronToField +const MeshField::LO tetNode2DofHolder[10] = {0, 1, 2, 3, 0, 2, 3, 1, 5, 4}; ``` -You must determine the correct rotation offset for your topology and DOF holders by -comparing the Omega\_h canonical vertex/edge ordering (see +You must determine the correct vertex and edge mapping for your topology and +DOF holders by comparing the Omega\_h canonical vertex/edge ordering (see figures above based on `Omega_h_simplex.hpp` and the `simplex_down_template` function) with the node ordering established by `getNodeParametricCoords()` in your shape struct. It is a good idea to add a unit test that evaluates the shape functions at the parametric coordinates of each node and checks that the value for that node is 1.0 and all others are 0.0. +### Jacobian Orientation + +@ref MeshField::FieldElement "FieldElement" returns **signed** Jacobian +determinants for 2D and 3D elements. The integrator +(`src/MeshField_Integrate.hpp`) uses the determinants as-is, and debug builds +assert that each one is non-negative. +The vertex ordering produced by your mapping must therefore follow the +right-hand rule so that the Jacobian determinant of a valid (un-inverted) +Omega\_h element is positive. +An incorrect vertex permutation results in a negative determinant and a +failed assertion in debug builds. + ### Interface Requirements | Member | Description | @@ -225,7 +249,7 @@ The steps within that branch are: 2. **Allocate storage** -- construct a `Controller` sized for each DOF-holding entity type. Each field in the controller needs entries for `numEntities * dofsPerHolder * numComp`. 3. **Create fields** -- call `MeshField::makeField(ctrl)` for each field index `N` to obtain typed slice objects. 4. **Construct the Accessor** -- pass the slice objects to your Accessor's constructor (e.g., `QuadraticAccessor{vtxField, edgeField}`). -5. **Return** `FieldWithController>{ctrl, ShapeField(meshInfo, accessor)}`. +5. **Return** `FieldWithController>{ctrl, ShapeField(meshInfo, accessor)}`. See the existing `order == 2 && dim == 2` branch (quadratic triangle) as a concrete reference: diff --git a/pumiTest/CMakeLists.txt b/pumiTest/CMakeLists.txt new file mode 100644 index 0000000..695a477 --- /dev/null +++ b/pumiTest/CMakeLists.txt @@ -0,0 +1,17 @@ +cmake_minimum_required(VERSION 3.16) +project(PumiTest LANGUAGES CXX C) + +find_package(MPI REQUIRED) + +find_package(SCOREC REQUIRED) + +add_executable(PumiTest main.cpp) + +target_link_libraries(PumiTest PRIVATE + MPI::MPI_CXX + SCOREC::apf + SCOREC::mds + SCOREC::gmi + SCOREC::pcu + SCOREC::lion +) diff --git a/pumiTest/main.cpp b/pumiTest/main.cpp new file mode 100644 index 0000000..9ae4bf2 --- /dev/null +++ b/pumiTest/main.cpp @@ -0,0 +1,56 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class CountIntegrator : public apf::Integrator { + protected: + double numEnt; + apf::MeshElement* currentMeshElement; + public: + double getCount() {return numEnt;} + void resetCount() { numEnt = 0; } + CountIntegrator(size_t order) : Integrator(order), numEnt(0) {}; + void inElement(apf::MeshElement* me) override { + currentMeshElement = me; + } + void atPoint(apf::Vector3 const& point, double weight, double jacobian) override { + apf::Vector3 globalCoord; + apf::mapLocalToGlobal(currentMeshElement, point, globalCoord); + numEnt += weight * globalCoord[0] * jacobian; + } +}; +int main(int argc, char ** argv) { + pcu::Init(&argc, &argv); + { + pcu::PCU PCUObj; + size_t order = atoi(argv[1]); + size_t size = atoi(argv[2]); + size_t numberOfElems = atoi(argv[3]); + size_t dim = atoi(argv[4]); + size_t runs = atoi(argv[5]); + gmi_register_mesh(); + gmi_register_null(); + apf::Mesh2* mesh = apf::makeMdsBox(numberOfElems, numberOfElems, dim == 3 ? numberOfElems : 0, size, size, dim == 3 ? size : 0, 1, &PCUObj); + CountIntegrator * countInt = new CountIntegrator(order); + for (int i = 0; i < runs; ++i) { + countInt->resetCount(); + auto start = std::chrono::high_resolution_clock::now(); + countInt->process(mesh); + auto end = std::chrono::high_resolution_clock::now(); + std::cout << "RESULT: " << pumi_git_sha() << "," << mesh->count(mesh->getDimension()) << "," << order << "," << countInt->getCount() << "," << std::chrono::duration(end - start).count() << std::endl; + } + delete countInt; + mesh->destroyNative(); + apf::destroyMesh(mesh); + } + pcu::Finalize(); + return 0; +} diff --git a/src/KokkosController.hpp b/src/KokkosController.hpp index 0d2deec..29d6ca3 100644 --- a/src/KokkosController.hpp +++ b/src/KokkosController.hpp @@ -62,7 +62,7 @@ template struct KokkosSliceWrapper { } }; -template +template class KokkosController { // type definitions @@ -182,7 +182,8 @@ class KokkosController { std::tuple...> values_; public: - typedef ExecutionSpace exe; + typedef ExeSpace exe; + using ExecutionSpace=exe; static const int MAX_RANK = 5; KokkosController() { diff --git a/src/MeshField.hpp b/src/MeshField.hpp index e9a38fc..a967837 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -53,6 +53,7 @@ struct SetCoordFieldFunctor const CoordinateMap map_; }; + template struct InterpolateElementFunctor { @@ -97,9 +98,7 @@ struct LinearTriangleToVertexField { const auto vtxDim = 0; const auto ignored = -1; const auto localVtxIdx = - (Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored) + - 2) % - 3; + Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored); const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); const MeshField::LO vtx = triVerts[(tri * triToVtxDegree) + localVtxIdx]; return {0, triCompIdx, vtx, MeshField::Vertex}; // node, comp, ent, topo @@ -127,10 +126,7 @@ struct LinearTetrahedronToVertexField { const auto tetDim = 3; const auto vtxDim = 0; const auto ignored = -1; - const auto localVtxIdx = - (Omega_h::simplex_down_template(tetDim, vtxDim, tetNodeIdx, ignored) + - 3) % - 4; + const auto localVtxIdx = Omega_h::simplex_down_template(tetDim, vtxDim, tetNodeIdx, ignored); const auto tetToVtxDegree = Omega_h::simplex_degree(tetDim, vtxDim); const MeshField::LO vtx = tetVerts[(tet * tetToVtxDegree) + localVtxIdx]; return {0, tetCompIdx, vtx, MeshField::Vertex}; // node, comp, ent, topo @@ -178,10 +174,8 @@ struct QuadraticTriangleToField { const auto triDim = 2; const auto vtxDim = 0; const auto ignored = -1; - const auto localVtxIdx = (Omega_h::simplex_down_template( - triDim, vtxDim, dofHolderIdx, ignored) + - 2) % - 3; + const auto localVtxIdx = Omega_h::simplex_down_template( + triDim, vtxDim, dofHolderIdx, ignored); const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx]; } else if (dofHolderTopo == MeshField::Edge) { @@ -191,7 +185,7 @@ struct QuadraticTriangleToField { // passing dofHolderIdx as Omega_h_simplex.hpp does not provide // a function that maps a triangle and edge index to a 'canonical' edge // index. This may need to be revisited... - osh_ent = triEdges[(tri * triToEdgeDegree) + (dofHolderIdx + 2) % 3]; + osh_ent = triEdges[(tri * triToEdgeDegree) + dofHolderIdx]; } else { assert(false); } @@ -222,7 +216,7 @@ struct QuadraticTetrahedronToField { operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx, MeshField::LO tet, MeshField::Mesh_Topology topo) const { assert(topo == MeshField::Tetrahedron); - const MeshField::LO tetNode2DofHolder[10] = {0, 1, 2, 3, 3, 4, 5, 0, 1, 2}; + const MeshField::LO tetNode2DofHolder[10] = {0, 1, 2, 3, 0, 2, 3, 1, 5, 4}; const MeshField::Mesh_Topology tetNode2DofHolderTopo[10] = { MeshField::Vertex, MeshField::Vertex, MeshField::Vertex, MeshField::Vertex, MeshField::Edge, MeshField::Edge, @@ -235,12 +229,8 @@ struct QuadraticTetrahedronToField { const auto tetDim = 3; const auto vtxDim = 0; const auto ignored = -1; - // cyclic rotation of the omegah vertex order to map to the meshfields order - // defined by the shape functions in MeshField_Shape.hpp - const auto localVtxIdx = (Omega_h::simplex_down_template( - tetDim, vtxDim, dofHolderIdx, ignored) + - 3) % - 4; + const auto localVtxIdx = Omega_h::simplex_down_template( + tetDim, vtxDim, dofHolderIdx, ignored); const auto tetToVtxDegree = Omega_h::simplex_degree(tetDim, vtxDim); osh_ent = tetVerts[(tet * tetToVtxDegree) + localVtxIdx]; } else if (dofHolderTopo == MeshField::Edge) { @@ -357,8 +347,8 @@ class OmegahMeshField { if constexpr (order == 1) { const auto meshDim = meshInfo.dim; const auto coords = mesh.coords(); - MeshField::parallel_for( - ExecutionSpace(), {0}, {meshInfo.numVtx}, + MeshField::parallel_for( + {0}, {meshInfo.numVtx}, SetCoordFieldFunctor(field, coords, meshDim, map), "setCoordField"); diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 1ab15d7..7e8fc55 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include // getLastValue #include @@ -48,9 +49,9 @@ template class has_static_rank { template class has_extent_method { private: template - static auto test(int) - -> decltype(std::declval().extent(std::declval()), - std::true_type()); + static auto + test(int) -> decltype(std::declval().extent(std::declval()), + std::true_type()); template static std::false_type test(...); public: @@ -308,8 +309,7 @@ struct FieldElement { // TODO use nested parallel for? auto Ji = Kokkos::subview(J, i, Kokkos::ALL(), Kokkos::ALL()); const auto cross = Ji(0, 0) * Ji(1, 1) - Ji(1, 0) * Ji(0, 1); - const auto magnitude = Kokkos::fabs(cross); - determinants(i) = magnitude; + determinants(i) = cross; }); return determinants; } @@ -326,6 +326,80 @@ struct FieldElement { return Kokkos::View("foo", J.extent(0)); } + Kokkos::View getJacobiansStaticLength(Kokkos::View localCoords) { + if (Debug) { + LO numErrors = 0; + Kokkos::parallel_reduce( + "checkCoords", localCoords.extent(0), + KOKKOS_LAMBDA(const int &ent, LO &lerrors) { + Real sum = 0; + LO isError = 0; + for (int i = 0; i < localCoords.extent(1); i++) { + if (localCoords(ent, i) < 0) + isError++; + sum += localCoords(ent, i); + } + if (Kokkos::fabs(sum - 1) > MachinePrecision) + isError++; + lerrors += isError; + }, + numErrors); + if (numErrors) { + fail("One or more of the parametric coordinates passed " + "to evaluate(...) were invalid\n"); + } + } + if (MeshEntDim != 1 && MeshEntDim != 2 && MeshEntDim != 3) { + fail("getJacobians only currently supports 1d, 2d, and 3d meshes. Input " + "mesh " + "has %zu dimensions.\n", + numMeshEnts); + } + + if constexpr (MeshEntDim == 1) { + const auto numPts = localCoords.extent(0); + Kokkos::View res("result", numPts * numMeshEnts, 1, 1); + auto jacobianFunc = KOKKOS_CLASS_LAMBDA(const int ent, const int pt) { + const auto val = getJacobian1d(ent); + res(pt, 0, 0) = val; + }; + MeshField::parallel_for({0, 0}, + {numMeshEnts, numPts}, jacobianFunc, "1dJacobian"); + return res; + } else if constexpr (MeshEntDim == 2 || MeshEntDim == 3) { + const auto numPts = localCoords.extent(0); + Kokkos::View res("result", numPts * numMeshEnts, MeshEntDim, + MeshEntDim); + Kokkos::deep_copy(res, 0.0); + auto jacobianFunc = KOKKOS_CLASS_LAMBDA(const int ent, const int pt) { + const auto vals = getNodeValues(ent); + Kokkos::Array xi; + for (size_t d = 0; d < MeshEntDim; d++) + xi[d] = localCoords(pt, d); + const auto grad = shapeFn.getLocalGradients(xi); + auto A = Kokkos::subview(res, ent * numPts + pt, Kokkos::ALL(), + Kokkos::ALL()); + Real localA[MeshEntDim][MeshEntDim] = {}; + for (size_t node = 0; node < ShapeType::numNodes; node++) { + for (size_t i = 0; i < MeshEntDim; ++i) { + for (size_t j = 0; j < MeshEntDim; ++j) { + localA[j][i] += vals[node * MeshEntDim + i] * + grad[node * MeshEntDim + j]; + } + } + } + for (size_t i = 0; i < MeshEntDim; ++i) { + for (size_t j = 0; j < MeshEntDim; ++j) { + A(i, j) = localA[i][j]; + } + } + }; + MeshField::parallel_for({0, 0}, + {numMeshEnts, numPts}, jacobianFunc, "2d3dJacobian"); + return res; + } + } + /** * @brief * Given an array of parametric coordinates 'localCoords', one per mesh @@ -404,48 +478,80 @@ struct FieldElement { // one matrix per point Kokkos::View res("result", numPts, MeshEntDim, MeshEntDim); Kokkos::deep_copy(res, 0.0); // initialize all entries to zero - - // fill the views of node coordinates and node gradients - Kokkos::View nodeCoords( - "nodeCoords", numPts); - Kokkos::View nodalGradients( - "nodalGradients", numPts); Kokkos::parallel_for( numMeshEnts, KOKKOS_CLASS_LAMBDA(const int ent) { const auto vals = getNodeValues(ent); assert(vals.size() == MeshEntDim * ShapeType::numNodes); for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) { + auto A = Kokkos::subview(res, pt, Kokkos::ALL(), Kokkos::ALL()); + Real localA[MeshEntDim][MeshEntDim] = {}; Kokkos::Array xi; for (size_t d = 0; d < MeshEntDim; d++) xi[d] = localCoords(pt, d); const auto grad = shapeFn.getLocalGradients(xi); for (size_t node = 0; node < ShapeType::numNodes; node++) { - for (size_t d = 0; d < MeshEntDim; d++) { - nodeCoords(pt, node, d) = vals[node * MeshEntDim + d]; - nodalGradients(pt, node, d) = grad[node * MeshEntDim + d]; + for (size_t i = 0; i < MeshEntDim; ++i) { + for (size_t j = 0; j < MeshEntDim; ++j) { + localA[j][i] += vals[node * MeshEntDim + i] * + grad[node * MeshEntDim + j]; + } } } - } - }); - - Kokkos::parallel_for( - numMeshEnts, KOKKOS_LAMBDA(const int ent) { - // TODO use nested parallel for? - for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) { - auto A = Kokkos::subview(res, pt, Kokkos::ALL(), Kokkos::ALL()); - for (size_t node = 0; node < ShapeType::numNodes; node++) { - auto a = - Kokkos::subview(nodalGradients, pt, node, Kokkos::ALL()); - auto b = Kokkos::subview(nodeCoords, pt, node, Kokkos::ALL()); - addTensorProduct(a, b, A); + for (size_t i = 0; i < MeshEntDim; ++i) { + for (size_t j = 0; j < MeshEntDim; ++j) { + A(i, j) = localA[i][j]; + } } } }); + return res; } } }; +template +Kokkos::View +evaluateStaticLength(FieldElement &fes, Kokkos::View localCoords) { + if (Debug) { + LO numErrors = 0; + Kokkos::parallel_reduce( + "checkCoords", localCoords.extent(0), + KOKKOS_LAMBDA(const int &ent, LO &lerrors) { + Real sum = 0; + LO isError = 0; + for (int i = 0; i < localCoords.extent(1); i++) { + if (localCoords(ent, i) < 0) + isError++; + sum += localCoords(ent, i); + } + if (Kokkos::fabs(sum - 1) > MachinePrecision) + isError++; + lerrors += isError; + }, + numErrors); + if (numErrors) { + fail("One or more of the parametric coordinates passed " + "to evaluate(...) were invalid\n"); + } + } + constexpr const auto numComponents = FieldElement::ValArray::size(); + const auto numPts = localCoords.extent(0); + Kokkos::View res("result", numPts * fes.numMeshEnts); + auto evaluateFunc = KOKKOS_LAMBDA(const int ent, const int pt) { + Kokkos::Array lc; + for (int i = 0; i < localCoords.extent(1); ++i) + lc[i] = localCoords(pt, i); + const auto val = fes.getValue(ent, lc); + for (int i = 0; i < numComponents; ++i) { + res(ent * numPts + pt, i) = val[i]; + } + }; + MeshField::parallel_for({0, 0}, {fes.numMeshEnts, numPts}, evaluateFunc, + "evaluate"); + return res; +} + /** * @brief * Given an array of parametric coordinates 'localCoords', one per mesh element, @@ -465,9 +571,9 @@ struct FieldElement { * @return Kokkos::View of evaluation results for all the mesh elements */ template -Kokkos::View -evaluate(FieldElement &fes, Kokkos::View localCoords, - Kokkos::View offsets) { +Kokkos::View evaluate( + FieldElement &fes, Kokkos::View localCoords, + Kokkos::View offsets) { if (Debug) { // check input parametric coords are positive and sum to one LO numErrors = 0; diff --git a/src/MeshField_For.hpp b/src/MeshField_For.hpp index 3b605ad..caeabd8 100644 --- a/src/MeshField_For.hpp +++ b/src/MeshField_For.hpp @@ -3,25 +3,44 @@ #include "MeshField_Utility.hpp" #include +#include "KokkosController.hpp" +#ifdef MESHFIELDS_ENABLE_CABANA +#include "MeshField_SimdFor.hpp" +#endif namespace MeshField { -template -void parallel_for(ExecutionSpace, const std::initializer_list &start, +template class Y> struct checkController { + static constexpr bool value = false; +}; + +template