Implement NonIsoparametricFieldElement - #99
Conversation
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
cwsmith
left a comment
cwsmith
left a comment
There was a problem hiding this comment.
Some comments are below. Mainly, I'm a bit concerned about the storage of the per-element coefficients into an single implicitly 'packed' view; I'd rather the contents be explicitly stored/accessed as it reduces the chance of indexing bugs and improves readability.
| #include "Omega_h_file.hpp" //move | ||
| #include "Omega_h_mesh.hpp" //move | ||
| #include "Omega_h_simplex.hpp" //move | ||
| #include <vector> |
| // These gradients are not currently used in the MeshField evaluation pipeline | ||
| // because the Reduced Quintic element uses a linear geometric mapping. The | ||
| // Jacobian is therefore computed from the linear geometry rather than from the | ||
| // Reduced Quintic shape functions. |
There was a problem hiding this comment.
Is there a use for this function? If not, I think we should remove this code. @jacobmerson thoughts?
| // xi = -b*λ0 + a*λ1 | ||
| // eta = c*λ2 | ||
| // λ0 + λ1 + λ2 = 1 |
There was a problem hiding this comment.
please remove/replace any non-ascii characters
| /** | ||
| * @brief Helper functions for reduced quintic coordinate transformations | ||
| */ | ||
| namespace ReducedQuinticHelpers { |
There was a problem hiding this comment.
I think moving these two functions into the ReducedQuinticTriangleShape struct as static member functions would improve the organization of the code. They are only be used within the struct (getValues and getLocalGradients) and within test/testReducedQuintic.cpp.
| * xi_local = a*λ1 - b*λ0 | ||
| * eta_local = c*λ2 |
There was a problem hiding this comment.
please remove/replace non-ascii characters
| * @param allDofs Array of all 6 DOFs for this vertex in physical coordinates | ||
| * @return The transformed DOF value in local coordinates | ||
| */ | ||
| template<typename Real> |
There was a problem hiding this comment.
why is the template needed?
| */ | ||
| template<typename Real> | ||
| KOKKOS_INLINE_FUNCTION | ||
| int solveLU_internal(int n, int nrhs, Real* A, int lda, Real* B, int ldb) { |
There was a problem hiding this comment.
I'd prefer we didn't use raw pointers as there is more chance for error etc. Can Kokkos Views be used instead so we can check their sizes etc.?
| * | ||
| * @param n Dimension of the matrix (20 for reduced quintic) | ||
| * @param nrhs Number of right-hand sides (18 for reduced quintic) | ||
| * @param A Input matrix (destroyed on output), row-major format |
There was a problem hiding this comment.
the 'destroyed on output' comment seems incorrect
| // Use fixed-size array for pivoting (n=20 maximum, which fits all our use cases) | ||
| // We use a stack array since this is device-compatible | ||
| int pivot[20]; |
There was a problem hiding this comment.
n should be asserted to be at most 20
| fail("mesh has no vertices\n"); | ||
| } | ||
| static constexpr size_t reducedQuinticNumComp = | ||
| 6; // value, dx, dy, dxx, dxy, dyy |
There was a problem hiding this comment.
It seems like this should be a constant encoded elsewhere (static member of ReducedQuinticTriangleShape?).
|
For "ReducedQuintic" we should consider if it makes sense to name Bell/Argyris which provides the more general family of constrained C1 elements. |
|
I suspect if we want to get performance, we need to look at the Kirby paper referenced here: #100 . My understanding (from the abstract) is that should allow us to not require a per-element style construction. He specifically mentions the Jardin paper. Note, the ReducedQuintic element we are working with is called the Bell element in that paper. |
| * @return 0 on success, positive value if singular | ||
| */ | ||
| KOKKOS_INLINE_FUNCTION | ||
| int solveLU_internal( |
There was a problem hiding this comment.
Rather than the current naming, we should put this into the detail namespace. Generally, I'd rather use a library for these sorts of functions because there are a lot of subtle points to getting these algorithms right. Since we already have a dependency on Kokkos, what about using KokkosKernels for this?
@cwsmith how do you feel about a kokkos kernels dependency here?
There was a problem hiding this comment.
On the fence. I like the idea of better performance. Is kokkos kernels already a dependency of pcms?
There was a problem hiding this comment.
Yes, KokkosKernels is a dependency for some of the field transfer methods (not all). M3DC1 also has an alternative "explicit" implementation which may not require the solver. We may want to look at the difference w.r.t. how they would perform on GPU.
4 participants
This is an alternative implementation of the reduced quintic element intended to replace and close #89 . The underlying implementation is the same, but the API has been extracted from the existing isoparametric field element API.
This approach may also potentially support other non-isoparametric elements in the future.