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 "integrator.hpp"
       5              : 
       6              : #include "fem/libceed/integrator.hpp"
       7              : 
       8              : namespace palace
       9              : {
      10              : 
      11              : namespace fem
      12              : {
      13              : 
      14        86071 : int DefaultIntegrationOrder::Get(const mfem::IsoparametricTransformation &T)
      15              : {
      16        86071 :   return 2 * p_trial + (q_order_jac ? T.OrderW() : 0) +
      17        86071 :          (T.GetFE()->Space() == mfem::FunctionSpace::Pk ? q_order_extra_pk
      18        86071 :                                                         : q_order_extra_qk);
      19              : }
      20              : 
      21         9870 : int DefaultIntegrationOrder::Get(const mfem::ElementTransformation &T)
      22              : {
      23         9870 :   const auto *T_iso = dynamic_cast<const mfem::IsoparametricTransformation *>(&T);
      24         9870 :   MFEM_VERIFY(
      25              :       T_iso,
      26              :       "Unexpected non-isoparametric element transformation to calculate quadrature order!");
      27         9870 :   return Get(*T_iso);
      28              : }
      29              : 
      30            0 : int DefaultIntegrationOrder::Get(const mfem::Mesh &mesh, mfem::Geometry::Type geom)
      31              : {
      32            0 :   MFEM_VERIFY(mesh.GetNodes(), "The mesh has no nodal FE space!");
      33            0 :   mfem::IsoparametricTransformation T;
      34            0 :   T.SetFE(mesh.GetNodalFESpace()->FEColl()->FiniteElementForGeometry(geom));
      35            0 :   return Get(T);
      36            0 : }
      37              : 
      38              : }  // namespace fem
      39              : 
      40         1254 : void DiscreteInterpolator::Assemble(Ceed ceed, CeedElemRestriction trial_restr,
      41              :                                     CeedElemRestriction test_restr, CeedBasis interp_basis,
      42              :                                     CeedOperator *op, CeedOperator *op_t)
      43              : {
      44              :   // Interpolators do not use an integration rule to map between the test and trial spaces.
      45         1254 :   ceed::AssembleCeedInterpolator(ceed, trial_restr, test_restr, interp_basis, op, op_t);
      46         1254 : }
      47              : 
      48          252 : void VectorFEBoundaryLFIntegrator::AssembleRHSElementVect(const mfem::FiniteElement &fe,
      49              :                                                           mfem::ElementTransformation &T,
      50              :                                                           mfem::Vector &elvect)
      51              : {
      52              :   const int dof = fe.GetDof();
      53              :   const int dim = fe.GetDim();
      54          252 :   const int q_order = fem::DefaultIntegrationOrder::Get(T);
      55          252 :   const mfem::IntegrationRule &ir = mfem::IntRules.Get(fe.GetGeomType(), q_order);
      56          252 :   f_hat.SetSize(dim);
      57          252 :   vshape.SetSize(dof, dim);
      58          252 :   elvect.SetSize(dof);
      59          252 :   elvect = 0.0;
      60              : 
      61         1008 :   for (int i = 0; i < ir.GetNPoints(); i++)
      62              :   {
      63              :     const mfem::IntegrationPoint &ip = ir.IntPoint(i);
      64              :     T.SetIntPoint(&ip);
      65          756 :     fe.CalcVShape(ip, vshape);
      66              : 
      67          756 :     Q.Eval(f_loc, T, ip);
      68          756 :     T.InverseJacobian().Mult(f_loc, f_hat);
      69         1512 :     f_hat *= ip.weight * T.Weight();
      70          756 :     vshape.AddMult(f_hat, elvect);
      71              :   }
      72          252 : }
      73              : 
      74            0 : void BoundaryLFIntegrator::AssembleRHSElementVect(const mfem::FiniteElement &fe,
      75              :                                                   mfem::ElementTransformation &T,
      76              :                                                   mfem::Vector &elvect)
      77              : {
      78              :   const int dof = fe.GetDof();
      79            0 :   const int q_order = fem::DefaultIntegrationOrder::Get(T);
      80            0 :   const mfem::IntegrationRule &ir = mfem::IntRules.Get(fe.GetGeomType(), q_order);
      81            0 :   shape.SetSize(dof);
      82            0 :   elvect.SetSize(dof);
      83            0 :   elvect = 0.0;
      84              : 
      85            0 :   for (int i = 0; i < ir.GetNPoints(); i++)
      86              :   {
      87              :     const mfem::IntegrationPoint &ip = ir.IntPoint(i);
      88              :     T.SetIntPoint(&ip);
      89            0 :     fe.CalcShape(ip, shape);
      90              : 
      91            0 :     double val = ip.weight * T.Weight() * Q.Eval(T, ip);
      92            0 :     elvect.Add(val, shape);
      93              :   }
      94            0 : }
      95              : 
      96              : }  // namespace palace