Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #include "coefficient.hpp"
       5              : 
       6              : #include <mfem.hpp>
       7              : #include "fem/libceed/ceed.hpp"
       8              : #include "models/materialoperator.hpp"
       9              : 
      10              : #include "fem/qfunctions/coeff/coeff_qf.h"
      11              : 
      12              : namespace palace::ceed
      13              : {
      14              : 
      15              : namespace
      16              : {
      17              : 
      18              : inline auto CoeffDim(int dim)
      19              : {
      20         9934 :   return dim * dim;
      21              : }
      22              : 
      23              : inline void MakeDiagonalCoefficient(int dim, CeedIntScalar *mat_coeff, CeedScalar a,
      24              :                                     CeedInt k)
      25              : {
      26              :   const int coeff_dim = CoeffDim(dim);
      27       189742 :   for (int i = 0; i < coeff_dim; i++)
      28              :   {
      29       161769 :     mat_coeff[coeff_dim * k + i].second = 0.0;
      30              :   }
      31        91600 :   for (int di = 0; di < dim; ++di)
      32              :   {
      33        63627 :     const int idx = di * (dim + 1);
      34        63627 :     mat_coeff[coeff_dim * k + idx].second = a;
      35              :   }
      36              : }
      37              : 
      38              : inline auto *AttrMat(CeedIntScalar *ctx)
      39              : {
      40              :   return ctx + 1;
      41              : }
      42              : 
      43              : inline auto *MatCoeff(CeedIntScalar *ctx)
      44              : {
      45       206974 :   const CeedInt num_attr = ctx[0].first;
      46       212519 :   return ctx + 2 + num_attr;
      47              : }
      48              : 
      49              : }  // namespace
      50              : 
      51        15479 : std::vector<CeedIntScalar> PopulateCoefficientContext(int dim,
      52              :                                                       const MaterialPropertyCoefficient *Q,
      53              :                                                       bool transpose, double a)
      54              : {
      55        15479 :   if (!Q)
      56              :   {
      57              :     // No attributes are stored in the map from attributes to material property coefficient,
      58              :     // indicating that all attributes map to the same identity coefficient.
      59         5545 :     std::vector<CeedIntScalar> ctx(2 + CoeffDim(dim), {0});
      60         5545 :     ctx[0].first = 0;
      61         5545 :     ctx[1].first = 1;
      62              :     MakeDiagonalCoefficient(dim, MatCoeff(ctx.data()), a, 0);
      63              :     return ctx;
      64              :   }
      65              : 
      66              :   // Material property coefficients might be empty if all attributes map to zero
      67              :   // coefficient.
      68              :   const auto &attr_mat = Q->GetAttributeToMaterial();
      69              :   const auto &mat_coeff = Q->GetMaterialProperties();
      70         9934 :   MFEM_VERIFY(attr_mat.Size() > 0, "Empty attributes for MaterialPropertyCoefficient!");
      71         9934 :   MFEM_VERIFY(attr_mat.Max() < mat_coeff.SizeK(),
      72              :               "Invalid attribute material property for MaterialPropertyCoefficient ("
      73              :                   << attr_mat.Max() << " vs. " << mat_coeff.SizeK() << ")!");
      74         9934 :   MFEM_VERIFY(mat_coeff.SizeI() == mat_coeff.SizeJ() &&
      75              :                   (mat_coeff.SizeI() == 1 || mat_coeff.SizeI() == dim),
      76              :               "Dimension mismatch for MaterialPropertyCoefficient and libCEED integrator!");
      77              : 
      78              :   // Map unassigned attributes to zero material property coefficient (the last material
      79              :   // property is reserved for zero).
      80              :   const int coeff_dim = CoeffDim(dim);
      81         9934 :   std::vector<CeedIntScalar> ctx(2 + attr_mat.Size() + coeff_dim * (mat_coeff.SizeK() + 1));
      82         9934 :   ctx[0].first = attr_mat.Size();
      83              :   const int zero_mat = mat_coeff.SizeK();
      84       144538 :   for (int i = 0; i < attr_mat.Size(); i++)
      85              :   {
      86       134604 :     const int k = attr_mat[i];
      87       269208 :     AttrMat(ctx.data())[i].first = (k < 0) ? zero_mat : k;
      88              :   }
      89              : 
      90              :   // Copy material properties.
      91         9934 :   ctx[1 + attr_mat.Size()].first = mat_coeff.SizeK() + 1;
      92        49436 :   for (int k = 0; k < mat_coeff.SizeK(); k++)
      93              :   {
      94        39502 :     if (mat_coeff.SizeI() == 1)
      95              :     {
      96              :       // Copy as diagonal matrix coefficient.
      97        22428 :       MakeDiagonalCoefficient(dim, MatCoeff(ctx.data()), a * mat_coeff(0, 0, k), k);
      98              :     }
      99              :     else
     100              :     {
     101        61824 :       for (int dj = 0; dj < dim; ++dj)
     102              :       {
     103       166056 :         for (int di = 0; di < dim; ++di)
     104              :         {
     105              :           // Column-major ordering.
     106       121306 :           const int idx = transpose ? (di * dim) + dj : (dj * dim) + di;
     107       121306 :           MatCoeff(ctx.data())[coeff_dim * k + idx].second = a * mat_coeff(di, dj, k);
     108              :         }
     109              :       }
     110              :     }
     111              :   }
     112        73174 :   for (int d = 0; d < coeff_dim; d++)
     113              :   {
     114        63240 :     MatCoeff(ctx.data())[coeff_dim * zero_mat + d].second = 0.0;
     115              :   }
     116              : 
     117              :   return ctx;
     118              : }
     119              : 
     120              : std::vector<CeedIntScalar>
     121           72 : PopulateCoefficientContext(int dim_mass, const MaterialPropertyCoefficient *Q_mass, int dim,
     122              :                            const MaterialPropertyCoefficient *Q, bool transpose_mass,
     123              :                            bool transpose, double a_mass, double a)
     124              : {
     125              :   // Mass coefficient comes first, then the other one for the QFunction.
     126           72 :   auto ctx_mass = PopulateCoefficientContext(dim_mass, Q_mass, transpose_mass, a_mass);
     127           72 :   auto ctx = PopulateCoefficientContext(dim, Q, transpose, a);
     128           72 :   ctx_mass.insert(ctx_mass.end(), ctx.begin(), ctx.end());
     129           72 :   return ctx_mass;
     130              : }
     131              : 
     132              : }  // namespace palace::ceed