Merge pull request #12734 from rppawlo/panzer-improve-mass-matrix-ass…

…embly

Panzer improve mass matrix assembly

packages/panzer/adapters-stk/src/Panzer_STK_LocalMeshUtilities.cpp

@@ -454,6 +454,8 @@ setupLocalMeshSidesetInfo(const panzer_stk::STK_Interface & mesh,
 Teuchos::RCP<panzer::LocalMeshInfo>
 generateLocalMeshInfo(const panzer_stk::STK_Interface & mesh)
 {
+  TEUCHOS_FUNC_TIME_MONITOR_DIFF("panzer_stk::generateLocalMeshInfo",GenerateLocalMeshInfo);
+
   using Teuchos::RCP;
   using Teuchos::rcp;
 
@@ -474,8 +476,6 @@ generateLocalMeshInfo(const panzer_stk::STK_Interface & mesh)
 
   Teuchos::RCP<const Teuchos::Comm<int> > comm = mesh.getComm();
 
-  TEUCHOS_FUNC_TIME_MONITOR_DIFF("panzer_stk::generateLocalMeshInfo",GenerateLocalMeshInfo);
-
   // This horrible line of code is required since the connection manager only takes rcps of a mesh
   RCP<const panzer_stk::STK_Interface> mesh_rcp = Teuchos::rcpFromRef(mesh);
   // We're allowed to do this since the connection manager only exists in this scope... even though it is also an RCP...

packages/panzer/adapters-stk/test/projection/projection.cpp

@@ -61,6 +61,7 @@
 #include "Panzer_Traits.hpp"
 #include "Panzer_Evaluator_WithBaseImpl.hpp"
 #include "Panzer_IntrepidBasisFactory.hpp"
+#include "Panzer_IntrepidOrientation.hpp"
 #include "Panzer_L2Projection.hpp"
 #include "Panzer_DOFManager.hpp"
 #include "Panzer_BlockedDOFManager.hpp"
@@ -152,7 +153,7 @@ TEUCHOS_UNIT_TEST(L2Projection, ToNodal)
   using LO = int;
   using GO = panzer::GlobalOrdinal;
   timer->start("ConnManager ctor");
-  const RCP<panzer::ConnManager> connManager = rcp(new panzer_stk::STKConnManager(mesh));
+  const RCP<panzer_stk::STKConnManager> connManager = rcp(new panzer_stk::STKConnManager(mesh));
   timer->stop("ConnManager ctor");
 
   // Set up bases for projections
@@ -244,6 +245,66 @@ TEUCHOS_UNIT_TEST(L2Projection, ToNodal)
   projectionFactory.setup(hgradBD,integrationDescriptor,comm,connManager,eBlockNames,worksetContainer);
   timer->stop("projectionFactory.setup()");
 
+
+
+
+
+  // Ignore the workset factory and set the basis directly
+  bool ignoreWorksetFactory = true;
+  std::map<std::string,Teuchos::RCP<panzer::BasisValues2<double>>> map_basis_values;
+  if (ignoreWorksetFactory) {
+
+    // Get the orientations for all element blocks
+    std::map<std::string,std::vector<Intrepid2::Orientation>> orientations;
+    panzer::buildIntrepidOrientations(eBlockNames,*connManager,orientations);
+
+    std::vector<shards::CellTopology> topologies;
+    connManager->getElementBlockTopologies(topologies);
+
+    for (size_t i=0; i <  eBlockNames.size(); ++i) {
+
+      const auto& eblock = eBlockNames[i];
+      const auto& topo = topologies[i];
+
+      // Get nodal coordinates for integration rules
+      std::vector<std::size_t> lids;
+      panzer::Intrepid2FieldPattern coord_fp(hgradBasis);
+      Kokkos::DynRankView<double,PHX::Device> nodal_coords_kokkos;
+      connManager->getDofCoords(eblock,coord_fp,lids,nodal_coords_kokkos);
+
+      // This is ugly. Need to fix connManager interface for layouts (fad vs non-fad).
+      PHX::MDField<double,panzer::Cell,panzer::BASIS,panzer::Dim> nodal_coords("nodal_coords","<Cell,BASIS,Dim>",
+                                                                               nodal_coords_kokkos.extent(0),
+                                                                               nodal_coords_kokkos.extent(1),
+                                                                               nodal_coords_kokkos.extent(2));
+      Kokkos::deep_copy(nodal_coords.get_view(),nodal_coords_kokkos);
+
+      // Build IV
+      Teuchos::RCP<shards::CellTopology> rcp_topo = Teuchos::rcp(new shards::CellTopology(topo.getCellTopologyData()));
+      Teuchos::RCP<panzer::IntegrationRule> ir = Teuchos::rcp(new panzer::IntegrationRule(integrationDescriptor,rcp_topo,static_cast<int>(nodal_coords_kokkos.extent(0))));
+      Teuchos::RCP<panzer::IntegrationValues2<double>> iv = Teuchos::rcp(new panzer::IntegrationValues2<double>);
+      iv->setup(ir,nodal_coords);
+
+      // Build BV
+      Teuchos::RCP<panzer::BasisValues2<double>> bv = Teuchos::rcp(new panzer::BasisValues2<double>);
+      Teuchos::RCP<panzer::BasisIRLayout> birl = Teuchos::rcp(new panzer::BasisIRLayout(hgradBD.getType(),hgradBD.getOrder(),*ir));
+      bv->setupUniform(birl,
+                       iv->getUniformCubaturePointsRef(true),
+                       iv->getJacobian(true),
+                       iv->getJacobianDeterminant(true),
+                       iv->getJacobianInverse(true));
+      bv->setOrientations(orientations[eblock]);
+      bv->setWeightedMeasure(iv->getWeightedMeasure(true));
+      bv->setCellNodeCoordinates(nodal_coords);
+      map_basis_values[eblock] = bv;
+    }
+    projectionFactory.useBasisValues(map_basis_values);
+  }
+
+
+
+
+
   TEST_ASSERT(nonnull(projectionFactory.getTargetGlobalIndexer()));
 
   // Build mass matrix

packages/panzer/disc-fe/src/Panzer_BasisDescriptor.cpp

@@ -134,6 +134,13 @@ getPointDescriptor() const
   return pd;
 }
 
+bool operator==(const panzer::BasisDescriptor& left,
+                const panzer::BasisDescriptor& right)
+{
+  return ( (left.getType() == right.getType()) &&
+           (left.getOrder() == right.getOrder()) );
+}
+
 } // end namespace panzer
 
 std::size_t

packages/panzer/disc-fe/src/Panzer_BasisDescriptor.hpp

@@ -66,8 +66,8 @@ class BasisDescriptor
 
   /** \brief Constructor for basis description
    *
-   * \param[in] basis_order Basis order (e.g. 1 could be piecewise linear)
-   * \param[in] basis_type Basis type (e.g. HGrad, HDiv, HCurl, ...)
+   * \param[in] basis_order Basis order as defined by Intrepid2 (e.g. 1 could be piecewise linear)
+   * \param[in] basis_type Basis type (a string: "HGrad", "HDiv", "HCurl", or "HVol")
    */
   BasisDescriptor(const int basis_order, const std::string & basis_type);
 
@@ -90,38 +90,35 @@ class BasisDescriptor
    */
   std::size_t getKey() const {return _key;}
 
-
   /** \brief Build a point descriptor that builds reference points for
    *  the DOF locations. This method throws if no points exist for this
    *  basis.
    */
   PointDescriptor getPointDescriptor() const;
 
 protected:
-
 
-  /// Basis type (HGrad, HDiv, HCurl,...)
+  /// Basis type (HGrad, HDiv, HCurl, HVol)
   std::string _basis_type;
 
   // Basis order (>0)
   int _basis_order;
 
   // Unique key associated with basis.
   std::size_t _key;
-
 };
 
-}
+bool operator==(const panzer::BasisDescriptor& left,
+                const panzer::BasisDescriptor& right);
 
+} // namespace panzer
 
 namespace std {
-
-template <>
-struct hash<panzer::BasisDescriptor>
-{
-  std::size_t operator()(const panzer::BasisDescriptor& desc) const;
-};
-
+  template <>
+  struct hash<panzer::BasisDescriptor>
+  {
+    std::size_t operator()(const panzer::BasisDescriptor& desc) const;
+  };
 }
 
 

packages/panzer/disc-fe/src/Panzer_BasisValues2.hpp

@@ -384,6 +384,9 @@ namespace panzer {
     hasUniformReferenceSpace() const
     {return is_uniform_;}
 
+    /// Return the basis descriptor
+    panzer::BasisDescriptor getBasisDescriptor() const;
+
     /// Get the extended dimensions used by sacado AD allocations
     const std::vector<PHX::index_size_type> &
     getExtendedDimensions() const