Initial commit for Exaconstit using libCEED through MFEM.

fem/bilininteg.hpp

@@ -1682,7 +1682,7 @@ private:
    const GeometricFactors *geom;  ///< Not owned
    int dim, ne, dofs1D, quad1D;
    Vector pa_data;
-   
+
 #ifdef MFEM_USE_CEED
    // CEED extension
    CeedData* ceedDataPtr;
@@ -1734,6 +1734,79 @@ public:
                                          const FiniteElement &test_fe);
 };
 
+/** Class for integrating the bilinear form a(u,v) := (Q grad u, grad v) where Q
+    can be a scalar or a matrix coefficient. */
+class MechanicsIntegrator: public BilinearFormIntegrator
+{
+protected:
+   Coefficient *Q;
+   MatrixCoefficient *MQ;
+
+private:
+   Vector vec, pointflux, shape;
+#ifndef MFEM_THREAD_SAFE
+   DenseMatrix dshape, dshapedxt, invdfdx, mq;
+   DenseMatrix te_dshape, te_dshapedxt;
+#endif
+
+   // PA extension
+   const DofToQuad *maps;         ///< Not owned
+   const GeometricFactors *geom;  ///< Not owned
+   int dim, ne, dofs1D, quad1D;
+   Vector pa_data;
+   
+#ifdef MFEM_USE_CEED
+   // CEED extension
+   CeedData* ceedDataPtr;
+   double* ktan_ptr;
+#endif
+
+public:
+   /// Construct a diffusion integrator with coefficient Q = 1
+   MechanicsIntegrator() { Q = NULL; MQ = NULL; maps = NULL; geom = NULL; }
+
+   /// Construct a diffusion integrator with a scalar coefficient q
+   MechanicsIntegrator(Coefficient &q)
+      : Q(&q) { MQ = NULL; maps = NULL; geom = NULL; }
+
+   /// Construct a diffusion integrator with a matrix coefficient q
+   MechanicsIntegrator(MatrixCoefficient &q)
+      : MQ(&q) { Q = NULL; maps = NULL; geom = NULL; }
+
+   /** Given a particular Finite Element
+       computes the element stiffness matrix elmat. */
+   virtual void AssembleElementMatrix(const FiniteElement &el,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat);
+   /** Given a trial and test Finite Element computes the element stiffness
+       matrix elmat. */
+   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                                       const FiniteElement &test_fe,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+
+   /// Perform the local action of the BilinearFormIntegrator
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &Tr,
+                                      const Vector &elfun, Vector &elvect);
+
+   virtual void ComputeElementFlux(const FiniteElement &el,
+                                   ElementTransformation &Trans,
+                                   Vector &u, const FiniteElement &fluxelem,
+                                   Vector &flux, int with_coef = 1);
+
+   virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
+                                    ElementTransformation &Trans,
+                                    Vector &flux, Vector *d_energy = NULL);
+
+   virtual void AssemblePA(const FiniteElementSpace&);
+
+   virtual void AddMultPA(const Vector&, Vector&) const;
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe);
+};
+
 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */
 class MassIntegrator: public BilinearFormIntegrator
 {

fem/libceed/mechanics.cpp

@@ -0,0 +1,128 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "ceed.hpp"
+#include "../../general/device.hpp"
+
+#include "mechanics.h"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+void CeedPAMechanicsAssemble(const FiniteElementSpace &fes,
+                             const mfem::IntegrationRule &irm,
+                             const double* ktan_ptr,
+                             CeedData& ceedData)
+{
+   Ceed ceed(internal::ceed);
+   mfem::Mesh *mesh = fes.GetMesh();
+   const int ir_order = irm.GetOrder();
+   const mfem::IntegrationRule &ir =
+      mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order);
+   CeedInt nqpts, nelem = mesh->GetNE(), dim = mesh->SpaceDimension();
+   mesh->EnsureNodes();
+   InitCeedTensorBasisAndRestriction(fes, ir, ceed, &ceedData.basis, &ceedData.restr);
+
+   const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace();
+   MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
+   InitCeedTensorBasisAndRestriction(*mesh_fes, ir, ceed, &ceedData.mesh_basis, &ceedData.mesh_restr);
+   CeedBasisGetNumQuadraturePoints(ceedData.basis, &nqpts);
+
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts * nelem, dim * dim * dim * dim, &ceedData.restr_i);
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts * nelem, 1, &ceedData.mesh_restr_i);
+
+   CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &ceedData.node_coords);
+   CeedVectorSetArray(ceedData.node_coords, CEED_MEM_HOST, CEED_USE_POINTER,
+                      mesh->GetNodes()->GetData());
+
+   CeedVectorCreate(ceed, nelem * nqpts * dim * dim * dim * dim, &ceedData.rho);
+
+   // Context data to be passed to the 'f_build_diff' Q-function.
+   ceedData.build_ctx.dim = mesh->Dimension();
+   ceedData.build_ctx.space_dim = mesh->SpaceDimension();
+
+   std::string mech_qf_file = GetCeedPath() + "/mechanics.h";
+   std::string mech_qf;
+
+   // Create the Q-function that builds the diff operator (i.e. computes its
+   // quadrature data) and set its context data.
+   mech_qf = mech_qf_file + ":f_build_mech";
+   CeedQFunctionCreateInterior(ceed, 1, f_build_mech,
+                               mech_qf.c_str(),
+                               &ceedData.build_qfunc);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "ktan", dim * dim * dim * dim,
+                          CEED_EVAL_NONE);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "dx", dim * dim, CEED_EVAL_GRAD);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "weights", 1, CEED_EVAL_WEIGHT);
+   CeedQFunctionAddOutput(ceedData.build_qfunc, "qd", dim * dim * dim * dim,
+                          CEED_EVAL_NONE);
+   CeedQFunctionSetContext(ceedData.build_qfunc, &ceedData.build_ctx,
+                           sizeof(ceedData.build_ctx));
+
+   // Create the operator that builds the quadrature data for the diff operator.
+   CeedOperatorCreate(ceed, ceedData.build_qfunc, NULL, NULL,
+                      &ceedData.build_oper);
+   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
+   if (mesh_fes->GetOrdering()==Ordering::byVDIM)
+   {
+      lmode = CEED_TRANSPOSE;
+   }
+   CeedVector ktan;
+   CeedVectorCreate(ceed, nelem * nqpts * dim * dim * dim * dim, &ktan);
+   CeedVectorSetArray(ktan, CEED_MEM_DEVICE, CEED_USE_POINTER, const_cast<double*>(ktan_ptr));
+   CeedOperatorSetField(ceedData.build_oper, "ktan", ceedData.restr_i,
+                        CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, ktan);
+   CeedOperatorSetField(ceedData.build_oper, "dx", ceedData.mesh_restr, lmode,
+                        ceedData.mesh_basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.build_oper, "weights", ceedData.mesh_restr_i,
+                        CEED_NOTRANSPOSE,
+                        ceedData.mesh_basis, CEED_VECTOR_NONE);
+   CeedOperatorSetField(ceedData.build_oper, "qd", ceedData.restr_i,
+                        CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
+
+   // Compute the quadrature data for the diff operator.
+   CeedOperatorApply(ceedData.build_oper, ceedData.node_coords, ceedData.rho,
+                     CEED_REQUEST_IMMEDIATE);
+
+   // Create the Q-function that defines the action of the diff operator.
+   mech_qf = mech_qf_file + ":f_apply_mech";
+   CeedQFunctionCreateInterior(ceed, 1, f_apply_mech,
+                               mech_qf.c_str(),
+                               &ceedData.apply_qfunc);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "u", dim * dim, CEED_EVAL_GRAD);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "qd", dim * dim * dim * dim,
+                         CEED_EVAL_NONE);
+   CeedQFunctionAddOutput(ceedData.apply_qfunc, "v", dim * dim, CEED_EVAL_GRAD);
+   CeedQFunctionSetContext(ceedData.apply_qfunc, &ceedData.build_ctx,
+                           sizeof(ceedData.build_ctx));
+
+   // Create the diff operator.
+   CeedOperatorCreate(ceed, ceedData.apply_qfunc, NULL, NULL, &ceedData.oper);
+   CeedOperatorSetField(ceedData.oper, "u", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.oper, "qd", ceedData.restr_i, CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, ceedData.rho);
+   CeedOperatorSetField(ceedData.oper, "v", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.u);
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.v);
+}
+
+#endif
+
+}

fem/libceed/mechanics.h

@@ -0,0 +1,181 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+/// A structure used to pass additional data to f_build_diff and f_apply_diff
+struct BuildContext { CeedInt dim, space_dim; CeedScalar coeff; };
+
+/// libCEED Q-function for building quadrature data for a diffusion operator with a constant coefficient
+CEED_QFUNCTION(f_build_mech)(void *ctx, CeedInt Q,
+                                   const CeedScalar *const *in, CeedScalar *const *out)
+{
+   BuildContext *bc = (BuildContext*)ctx;
+   // in[0] is ktan
+   // in[1] is Jacobians with shape [dim, nc=dim, Q]
+   // in[2] is quadrature weights, size (Q)
+   //
+   // At every quadrature point, compute qw/det(J).adj(J).adj(J)^T and store
+   // the symmetric part of the result.
+   const CeedScalar *ktan = in[0], *J = in[1], *qw = in[2];
+   CeedScalar *qd = out[0];
+   switch (bc->dim + 10 * bc->space_dim)
+   {
+      case 11:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            qd[i] = ktan[i] * qw[i] / J[i];
+         }
+         break;
+      case 22:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            // J: 0 2   qd: 0 1   adj(J):  J22 -J12
+            //    1 3       1 2           -J21  J11
+            // const CeedScalar J11 = J[i + Q * 0];
+            // const CeedScalar J21 = J[i + Q * 1];
+            // const CeedScalar J12 = J[i + Q * 2];
+            // const CeedScalar J22 = J[i + Q * 3];
+            // const CeedScalar w = qw[i] / (J11 * J22 - J21 * J12);
+            // qd[i + Q * 0] =   coeff * w * (J12 * J12 + J22 * J22);
+            // qd[i + Q * 1] = - coeff * w * (J11 * J12 + J21 * J22);
+            // qd[i + Q * 2] =   coeff * w * (J11 * J11 + J21 * J21);
+            //TODO
+         }
+         break;
+      case 33:
+         for (CeedInt q = 0; q < Q; q++)
+         {
+            // J: 0 3 6   qd: 0 1 2
+            //    1 4 7       1 3 4
+            //    2 5 8       2 4 5
+            const CeedScalar J11 = J[q + Q * 0];
+            const CeedScalar J21 = J[q + Q * 1];
+            const CeedScalar J31 = J[q + Q * 2];
+            const CeedScalar J12 = J[q + Q * 3];
+            const CeedScalar J22 = J[q + Q * 4];
+            const CeedScalar J32 = J[q + Q * 5];
+            const CeedScalar J13 = J[q + Q * 6];
+            const CeedScalar J23 = J[q + Q * 7];
+            const CeedScalar J33 = J[q + Q * 8];
+            const CeedScalar A11 = J22 * J33 - J23 * J32;
+            const CeedScalar A12 = J13 * J32 - J12 * J33;
+            const CeedScalar A13 = J12 * J23 - J13 * J22;
+            const CeedScalar A21 = J23 * J31 - J21 * J33;
+            const CeedScalar A22 = J11 * J33 - J13 * J31;
+            const CeedScalar A23 = J13 * J21 - J11 * J23;
+            const CeedScalar A31 = J21 * J32 - J22 * J31;
+            const CeedScalar A32 = J12 * J31 - J11 * J32;
+            const CeedScalar A33 = J11 * J22 - J12 * J21;
+            const CeedScalar w = qw[q] / (J11 * A11 + J21 * A12 + J31 * A13);
+            // Load ktan
+            CeedScalar K[3][3][3][3];
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  for (int l = 0; l < 3; ++l) {
+                     for (int m = 0; m < 3; ++m) {
+                        K[j][k][l][m] = ktan [q + (j + k*3 + l*3*3 + m*3*3*3) * Q];
+                     }
+                  }
+               }
+            }
+            // ktan*J^-1
+            CeedScalar tmp[3][3][3][3];
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  for (int l = 0; l < 3; ++l) {
+                     tmp[j][k][l][0] = K[j][k][l][0] * A11 + K[j][k][l][1] * A21 + K[j][k][l][2] * A31;
+                     tmp[j][k][l][1] = K[j][k][l][0] * A12 + K[j][k][l][1] * A22 + K[j][k][l][2] * A32;
+                     tmp[j][k][l][2] = K[j][k][l][0] * A13 + K[j][k][l][1] * A23 + K[j][k][l][2] * A33;
+                  }
+               }
+            }
+            // J^-T*ktan*J^-1
+            for (int k = 0; k < 3; ++k) {
+               for (int l = 0; l < 3; ++l) {
+                  for (int n = 0; n < 3; ++n) {
+                     qd[q + (0 + k*3 + l*3*3 + n*3*3*3)*Q] = w * (A11 * tmp[0][k][l][n] + A21 * tmp[1][k][l][n] + A31 * tmp[2][k][l][n]);
+                     qd[q + (1 + k*3 + l*3*3 + n*3*3*3)*Q] = w * (A12 * tmp[0][k][l][n] + A22 * tmp[1][k][l][n] + A32 * tmp[2][k][l][n]);
+                     qd[q + (2 + k*3 + l*3*3 + n*3*3*3)*Q] = w * (A13 * tmp[0][k][l][n] + A23 * tmp[1][k][l][n] + A33 * tmp[2][k][l][n]);
+                  }
+               }
+            }
+         }
+         break;
+   }
+   return 0;
+}
+
+/// libCEED Q-function for applying a diff operator
+CEED_QFUNCTION(f_apply_mech)(void *ctx, CeedInt Q,
+                             const CeedScalar *const *in, CeedScalar *const *out)
+{
+   BuildContext *bc = (BuildContext *)ctx;
+   // in[0], out[0] have shape [dim, nc=1, Q]
+   const CeedScalar *ug = in[0], *qd = in[1];
+   CeedScalar *vg = out[0];
+   switch (bc->dim)
+   {
+      case 1:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            vg[i] = ug[i] * qd[i];
+         }
+         break;
+      case 2:
+         // for (CeedInt i = 0; i < Q; i++)
+         // {
+         //    const CeedScalar ug0 = ug[i + Q * 0];
+         //    const CeedScalar ug1 = ug[i + Q * 1];
+         //    vg[i + Q * 0] = qd[i + Q * 0] * ug0 + qd[i + Q * 1] * ug1;
+         //    vg[i + Q * 1] = qd[i + Q * 1] * ug0 + qd[i + Q * 2] * ug1;
+         // }
+         break;
+      case 3:
+         for (CeedInt q = 0; q < Q; q++)
+         {
+            // Read spatial derivatives of u components
+            const CeedScalar uJ[3][3]        = {{ug[q+(0+0*3)*Q],
+                                                 ug[q+(0+1*3)*Q],
+                                                 ug[q+(0+2*3)*Q]},
+                                                {ug[q+(1+0*3)*Q],
+                                                 ug[q+(1+1*3)*Q],
+                                                 ug[q+(1+2*3)*Q]},
+                                                {ug[q+(2+0*3)*Q],
+                                                 ug[q+(2+1*3)*Q],
+                                                 ug[q+(2+2*3)*Q]}
+                                                };
+            // Load quadrature data
+            CeedScalar K[3][3][3][3];
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  for (int l = 0; l < 3; ++l) {
+                     for (int m = 0; m < 3; ++m) {
+                        K[j][k][l][m] = qd [q + (j + k*3 + l*3*3 + m*3*3*3) * Q];
+                     }
+                  }
+               }
+            }
+            // double contraction
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  vg[q + (j+k*3)*Q] = 0.0;
+                  for (int l = 0; l < 3; ++l) {
+                     for (int m = 0; m < 3; ++m) {
+                        vg[q + (j+k*3)*Q] += K[j][k][l][m] * uJ[m][l];
+                     }
+                  }
+               }
+            }
+         }
+         break;
+   }
+   return 0;
+}
+

fem/libceed/mechanics.hpp

@@ -0,0 +1,33 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_LIBCEED_MECH_HPP
+#define MFEM_LIBCEED_MECH_HPP
+
+#include "ceed.hpp"
+#include "../fespace.hpp"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+/// Initialize a Mechanics Integrator using libCEED
+void CeedPAMechanicsAssemble(const FiniteElementSpace &fes,
+                             const mfem::IntegrationRule &irm,
+                             const double* ktan_ptr,
+                             CeedData& ceedData);
+
+#endif
+
+}
+
+#endif // MFEM_LIBCEED_NECH_HPP