Merge pull request #1076 from sandialabs/bartgol/extruded-mesh-work

src/CMakeLists.txt

@@ -232,6 +232,9 @@ list (APPEND SOURCES
   disc/Albany_MeshSpecs.cpp
   disc/Albany_DOFManager.cpp
   disc/Albany_DiscretizationUtils.cpp
+  disc/Albany_ExtrudedMesh.cpp
+  disc/Albany_ExtrudedConnManager.cpp
+  disc/Albany_ExtrudedDiscretization.cpp
   )
 list (APPEND HEADERS
   disc/Albany_DiscretizationUtils.hpp
@@ -243,6 +246,9 @@ list (APPEND HEADERS
   disc/Albany_MeshSpecs.hpp
   disc/Albany_ConnManager.hpp
   disc/Albany_DOFManager.hpp
+  disc/Albany_ExtrudedMesh.hpp
+  disc/Albany_ExtrudedConnManager.hpp
+  disc/Albany_ExtrudedDiscretization.hpp
   )
 
 #stk

src/disc/Albany_AbstractDiscretization.hpp

@@ -13,6 +13,7 @@
 #include "Albany_DiscretizationUtils.hpp"
 #include "Albany_StateInfoStruct.hpp"
 #include "Albany_DOFManager.hpp"
+#include "Albany_ThyraCrsMatrixFactory.hpp"
 
 #include "Albany_ThyraTypes.hpp"
 #include "Albany_GlobalLocalIndexer.hpp"
@@ -134,8 +135,10 @@ class AbstractDiscretization
   }
 
   //! Create a Jacobian operator
-  virtual Teuchos::RCP<Thyra_LinearOp>
-  createJacobianOp() const = 0;
+  Teuchos::RCP<Thyra_LinearOp> createJacobianOp() const
+  {
+    return m_jac_factory->createOp();
+  }
 
   //! Get Node set lists
   virtual const NodeSetList&
@@ -243,8 +246,9 @@ class AbstractDiscretization
   }
 
   //! Get nodal parameters state info struct
-  virtual const StateInfoStruct&
-  getNodalParameterSIS() const = 0;
+  const StateInfoStruct& getNodalParameterSIS() const {
+    return getMeshStruct()->get_field_accessor()->getNodalParameterSIS();
+  }
 
   //! Retrieve Vector (length num worksets) of element block names
   const WorksetArray<std::string>&
@@ -388,6 +392,9 @@ class AbstractDiscretization
 protected:
   strmap_t<Teuchos::RCP<AbstractDiscretization>> sideSetDiscretizations;
 
+  //! Jacobian matrix operator factory
+  Teuchos::RCP<ThyraCrsMatrixFactory> m_jac_factory;
+
   // One dof mgr per dof per part
   // Notice that the dof mgr on a side is not the restriction
   // of the volume dof mgr to that side, since local ids are different.

src/disc/Albany_AbstractMeshStruct.hpp

@@ -44,6 +44,11 @@ struct AbstractMeshStruct {
 
     bool isBulkDataSet () const { return m_bulk_data_set; }
 
+    virtual LO get_num_local_elements () const = 0;
+    virtual LO get_num_local_nodes () const = 0;
+    virtual GO get_max_node_gid () const = 0;
+    virtual GO get_max_elem_gid () const = 0;
+
     virtual Teuchos::RCP<AbstractMeshFieldAccessor> get_field_accessor() const = 0;
 
     Teuchos::RCP<LayeredMeshNumbering<GO> > global_cell_layers_data;

src/disc/Albany_ConnManager.hpp

@@ -9,7 +9,8 @@
 
 #include "Albany_ScalarOrdinalTypes.hpp"
 
-#include "Panzer_ConnManager.hpp"
+#include <Panzer_ConnManager.hpp>
+#include <Panzer_FieldPattern.hpp>
 
 #include <vector>
 
@@ -100,10 +101,29 @@ class ConnManager : public panzer::ConnManager {
 
     return topologies[0];
   }
+
+  // Overload, not shadow
+  using panzer::ConnManager::buildConnectivity;
+  void buildConnectivity(const Teuchos::RCP<const panzer::FieldPattern> &fp) {
+    if (is_connectivity_built()) {
+      TEUCHOS_TEST_FOR_EXCEPTION (not m_fp->equals(*fp), std::runtime_error,
+          "Error! Rebuilding conn mgr with a different field pattern!\n"
+          "Old FP\n" << *m_fp << "\n"
+          "New FP\n" << *fp << "\n");
+      return;
+    }
+
+    this->buildConnectivity(*fp);
+
+    // Store copy of input pattern for later checks
+    m_fp = fp;
+  }
+
+  bool is_connectivity_built () const { return not m_fp.is_null(); }
 protected:
   std::vector<std::string> m_elem_blocks_names;
 
-  bool m_is_connectivity_built = false;
+  Teuchos::RCP<const panzer::FieldPattern> m_fp;
 };
 
 } // namespace Albany

src/disc/Albany_DOFManager.cpp

@@ -368,7 +368,7 @@ albanyBuildGlobalUnknowns ()
   // IMPORTANT! Do not use the GeometricAggFieldPattern, since you
   // *need* to count the same geo node multiple times if there are
   // multiple fields that need it.
-  connMngr_->buildConnectivity(*fa_fps_.back());
+  m_conn_mgr->buildConnectivity(fa_fps_.back());
 
   // We take a set as well, since std::find is on avg O(1) for unordered_set, vs O(N) in an array
   auto add_if_not_there = [](std::vector<GO>& v, std::unordered_set<GO>& s, const GO gid) {

src/disc/Albany_DOFManager.hpp

@@ -91,7 +91,7 @@ class DOFManager : public panzer::DOFManager {
     return m_conn_mgr->get_topology();
   }
 
-  Teuchos::RCP<const ConnManager> getAlbanyConnManager() const {
+  Teuchos::RCP<ConnManager> getAlbanyConnManager() const {
     return m_conn_mgr;
   }
 
@@ -146,7 +146,7 @@ class DOFManager : public panzer::DOFManager {
   // NOTE: this 
   vec4int       m_side_closure_orderd_as_side;
 
-  Teuchos::RCP<const ConnManager>           m_conn_mgr;
+  Teuchos::RCP<ConnManager>           m_conn_mgr;
 
   std::string m_part_name;
 

src/disc/Albany_DiscretizationFactory.cpp

@@ -8,6 +8,7 @@
 #include "Teuchos_TestForException.hpp"
 #include "Albany_DiscretizationFactory.hpp"
 
+#include "Albany_ExtrudedDiscretization.hpp"
 #include "Albany_STKDiscretization.hpp"
 // #include "Albany_BlockedSTKDiscretization.hpp"
 #include "Albany_TmplSTKMeshStruct.hpp"
@@ -18,6 +19,7 @@
 #include "Albany_AsciiSTKMesh2D.hpp"
 #include "Albany_GmshSTKMeshStruct.hpp"
 #include "Albany_ExtrudedSTKMeshStruct.hpp"
+#include "Albany_ExtrudedMesh.hpp"
 #include "Albany_Utils.hpp" // For CalculateNumberParams
 
 #ifdef ALBANY_OMEGAH
@@ -118,6 +120,30 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP<Teuchos::ParameterList> dis
     else if (method == "Extruded") {
         Teuchos::RCP<AbstractMeshStruct> basalMesh;
 
+        // Get basal_params
+        auto ss_disc_params = Teuchos::sublist(disc_params,"Side Set Discretizations");
+        auto basal_params = Teuchos::sublist(ss_disc_params,"basalside");
+        if (!basal_params->isParameter("Number Of Time Derivatives")) {
+          basal_params->set("Number Of Time Derivatives",disc_params->get<int>("Number Of Time Derivatives"));
+        }
+
+        // Set basal workset size
+        int extruded_ws_size = disc_params->get("Workset Size", -1);
+        if (extruded_ws_size == -1) {
+          basal_params->set("Workset Size", -1);
+        } else if (!basal_params->isParameter("Workset Size")) {
+          // Compute basal workset size based on extruded workset size
+          int basal_ws_size = extruded_ws_size / disc_params->get<int>("NumLayers");
+          basal_ws_size = std::max(basal_ws_size,1); //makes sure is at least 1.
+          basal_params->set("Workset Size", basal_ws_size);
+        }
+
+        basalMesh = createMeshStruct(basal_params, comm, numParams);
+        mesh = Teuchos::rcp(new ExtrudedMesh(basalMesh, disc_params, comm));
+    }
+    else if (method == "STKExtruded") {
+        Teuchos::RCP<AbstractMeshStruct> basalMesh;
+
         // Get basal_params
         Teuchos::RCP<Teuchos::ParameterList> basal_params;
         if (disc_params->isSublist("Side Set Discretizations") && disc_params->sublist("Side Set Discretizations").isSublist("basalside")) {
@@ -154,7 +180,7 @@ DiscretizationFactory::createMeshStruct(Teuchos::RCP<Teuchos::ParameterList> dis
                   "!" << std::endl << "Supplied parameter list is " << std::endl << *disc_params <<
                   "\nValid Methods are: STK1D, STK2D, STK3D, STK3DPoint, Ioss," <<
                   " Exodus, Ascii," <<
-                  " Ascii2D, Extruded" << std::endl);
+                  " Ascii2D, STKExtruded, Extruded" << std::endl);
   }
 
   if (disc_params->isSublist ("Side Set Discretizations")) {
@@ -304,7 +330,13 @@ createDiscretizationFromMeshStruct (const Teuchos::RCP<AbstractMeshStruct>& mesh
       rigidBodyModes->setPiroPL(piroParams);
 
   Teuchos::RCP<AbstractDiscretization> disc;
-  if (mesh->meshLibName()=="STK") {
+  if (mesh->meshSpecs[0]->mesh_type==MeshType::Extruded)
+  {
+    auto ext_mesh = Teuchos::rcp_dynamic_cast<ExtrudedMesh>(mesh);
+    auto basal_mesh = ext_mesh->basal_mesh();
+    auto basal_disc = createDiscretizationFromMeshStruct(basal_mesh,neq,{},rigidBodyModes);
+    disc = Teuchos::rcp(new ExtrudedDiscretization (discParams,neq,ext_mesh,basal_disc,comm,rigidBodyModes, sideSetEquations));
+  } else if (mesh->meshLibName()=="STK") {
     auto ms = Teuchos::rcp_dynamic_cast<AbstractSTKMeshStruct>(mesh);
     disc = Teuchos::rcp(new STKDiscretization(discParams, neq, ms, comm, rigidBodyModes, sideSetEquations));
 #ifdef ALBANY_OMEGAH