Albany: a Trilinos-based PDE code: Albany_MpasSTKMeshStruct.cpp Source File

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include <iostream>
00008 
00009 #include "Albany_MpasSTKMeshStruct.hpp"
00010 #include "Teuchos_VerboseObject.hpp"
00011 
00012 #include <Shards_BasicTopologies.hpp>
00013 
00014 #include <stk_mesh/base/Entity.hpp>
00015 #include <stk_mesh/base/GetEntities.hpp>
00016 #include <stk_mesh/base/GetBuckets.hpp>
00017 #include <stk_mesh/base/FieldData.hpp>
00018 #include <stk_mesh/base/Selector.hpp>
00019 #include <stk_io/IossBridge.hpp>
00020 #include <Ioss_SubSystem.h>
00021 
00022 //#include <stk_mesh/fem/FEMHelpers.hpp>
00023 #include <boost/algorithm/string/predicate.hpp>
00024 
00025 #include "Albany_Utils.hpp"
00026 
00027 Albany::MpasSTKMeshStruct::MpasSTKMeshStruct(const Teuchos::RCP<Teuchos::ParameterList>& params,
00028                                              const Teuchos::RCP<const Epetra_Comm>& comm,
00029                                              const std::vector<int>& indexToTriangleID, const std::vector<int>& verticesOnTria, int nGlobalTriangles) :
00030   GenericSTKMeshStruct(params,Teuchos::null, 2),
00031   out(Teuchos::VerboseObjectBase::getDefaultOStream()),
00032   periodic(false),
00033   NumEles(indexToTriangleID.size()),
00034   hasRestartSol(false),
00035   restartTime(0.)
00036 {
00037   elem_map = Teuchos::rcp(new Epetra_Map(nGlobalTriangles, indexToTriangleID.size(), &indexToTriangleID[0], 0, *comm)); // Distribute the elems equally
00038   
00039   params->validateParameters(*getValidDiscretizationParameters(),0);
00040 
00041 
00042   std::string ebn="Element Block 0";
00043   partVec[0] = & metaData->declare_part(ebn, metaData->element_rank() );
00044   ebNameToIndex[ebn] = 0;
00045 
00046 #ifdef ALBANY_SEACAS
00047   stk::io::put_io_part_attribute(*partVec[0]);
00048 #endif
00049 
00050 
00051   std::vector<std::string> nsNames;
00052   std::string nsn="Lateral";
00053   nsNames.push_back(nsn);
00054   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00055 #ifdef ALBANY_SEACAS
00056     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00057 #endif
00058   nsn="Internal";
00059   nsNames.push_back(nsn);
00060   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00061 #ifdef ALBANY_SEACAS
00062     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00063 #endif
00064 
00065 
00066   std::vector<std::string> ssNames;
00067   std::string ssn="LateralSide";
00068   ssNames.push_back(ssn);
00069     ssPartVec[ssn] = & metaData->declare_part(ssn, metaData->side_rank() );
00070 #ifdef ALBANY_SEACAS
00071     stk::io::put_io_part_attribute(*ssPartVec[ssn]);
00072 #endif
00073 
00074   stk::mesh::fem::set_cell_topology<shards::Triangle<3> >(*partVec[0]);
00075   stk::mesh::fem::set_cell_topology<shards::Line<2> >(*ssPartVec[ssn]);
00076 
00077   numDim = 2;
00078   int cub = params->get("Cubature Degree",3);
00079   int worksetSizeMax = params->get("Workset Size",50);
00080   int worksetSize = this->computeWorksetSize(worksetSizeMax, elem_map->NumMyElements());
00081 
00082   const CellTopologyData& ctd = *metaData->get_cell_topology(*partVec[0]).getCellTopologyData();
00083 
00084   this->meshSpecs[0] = Teuchos::rcp(new Albany::MeshSpecsStruct(ctd, numDim, cub,
00085                              nsNames, ssNames, worksetSize, partVec[0]->name(),
00086                              ebNameToIndex, this->interleavedOrdering));
00087 
00088 
00089 }
00090 
00091 //Wedge
00092 Albany::MpasSTKMeshStruct::MpasSTKMeshStruct(const Teuchos::RCP<Teuchos::ParameterList>& params,
00093                                              const Teuchos::RCP<const Epetra_Comm>& comm,
00094                                              const std::vector<int>& indexToTriangleID, const std::vector<int>& verticesOnTria, int nGlobalTriangles, int numLayers, int Ordering) :
00095   GenericSTKMeshStruct(params,Teuchos::null,3),
00096   out(Teuchos::VerboseObjectBase::getDefaultOStream()),
00097   periodic(false),
00098   NumEles(indexToTriangleID.size()),
00099   hasRestartSol(false),
00100   restartTime(0.)
00101 {
00102   std::vector<int> indexToPrismID(indexToTriangleID.size()*numLayers);
00103 
00104   //Int ElemColumnShift = (ordering == ColumnWise) ? 1 : indexToTriangleID.size();
00105   int elemColumnShift = (Ordering == 1) ? 1 : nGlobalTriangles;
00106   int lElemColumnShift = (Ordering == 1) ? 1 : indexToTriangleID.size();
00107   int elemLayerShift = (Ordering == 0) ? 1 : numLayers;
00108 
00109   for(int il=0; il< numLayers; il++)
00110   {
00111     int shift = il*elemColumnShift;
00112     int lShift = il*lElemColumnShift;
00113     for(int j=0; j< indexToTriangleID.size(); j++)
00114     {
00115       int lid = lShift + j*elemLayerShift;
00116       indexToPrismID[lid] = shift+elemLayerShift * indexToTriangleID[j];
00117     }
00118   }
00119 
00120   elem_map = Teuchos::rcp(new Epetra_Map(nGlobalTriangles*numLayers, indexToPrismID.size(), &indexToPrismID[0], 0, *comm)); // Distribute the elems equally
00121 
00122   params->validateParameters(*getValidDiscretizationParameters(),0);
00123 
00124 
00125   std::string ebn="Element Block 0";
00126   partVec[0] = & metaData->declare_part(ebn, metaData->element_rank() );
00127   ebNameToIndex[ebn] = 0;
00128 
00129 #ifdef ALBANY_SEACAS
00130   stk::io::put_io_part_attribute(*partVec[0]);
00131 #endif
00132 
00133 
00134   std::vector<std::string> nsNames;
00135   std::string nsn="Lateral";
00136   nsNames.push_back(nsn);
00137   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00138 #ifdef ALBANY_SEACAS
00139     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00140 #endif
00141   nsn="Internal";
00142   nsNames.push_back(nsn);
00143   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00144 #ifdef ALBANY_SEACAS
00145     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00146 #endif
00147   nsn="Bottom";
00148   nsNames.push_back(nsn);
00149   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00150 #ifdef ALBANY_SEACAS
00151   stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00152 #endif
00153 
00154 
00155   std::vector<std::string> ssNames;
00156   std::string ssnLat="lateralside";
00157   std::string ssnBottom="basalside";
00158   std::string ssnTop="upperside";
00159 
00160   ssNames.push_back(ssnLat);
00161   ssNames.push_back(ssnBottom);
00162   ssNames.push_back(ssnTop);
00163   ssPartVec[ssnLat] = & metaData->declare_part(ssnLat, metaData->side_rank() );
00164   ssPartVec[ssnBottom] = & metaData->declare_part(ssnBottom, metaData->side_rank() );
00165   ssPartVec[ssnTop] = & metaData->declare_part(ssnTop, metaData->side_rank() );
00166 #ifdef ALBANY_SEACAS
00167     stk::io::put_io_part_attribute(*ssPartVec[ssnLat]);
00168     stk::io::put_io_part_attribute(*ssPartVec[ssnBottom]);
00169     stk::io::put_io_part_attribute(*ssPartVec[ssnTop]);
00170 #endif
00171 
00172   stk::mesh::fem::set_cell_topology<shards::Wedge<6> >(*partVec[0]);
00173   stk::mesh::fem::set_cell_topology<shards::Triangle<3> >(*ssPartVec[ssnBottom]);
00174   stk::mesh::fem::set_cell_topology<shards::Triangle<3> >(*ssPartVec[ssnTop]);
00175   stk::mesh::fem::set_cell_topology<shards::Quadrilateral<4> >(*ssPartVec[ssnLat]);
00176 
00177   numDim = 3;
00178   int cub = params->get("Cubature Degree",3);
00179   int worksetSizeMax = params->get("Workset Size",50);
00180   int worksetSize = this->computeWorksetSize(worksetSizeMax, elem_map->NumMyElements());
00181 
00182   const CellTopologyData& ctd = *metaData->get_cell_topology(*partVec[0]).getCellTopologyData();
00183 
00184   this->meshSpecs[0] = Teuchos::rcp(new Albany::MeshSpecsStruct(ctd, numDim, cub,
00185                              nsNames, ssNames, worksetSize, partVec[0]->name(),
00186                              ebNameToIndex, this->interleavedOrdering));
00187 
00188 
00189 }
00190 
00191 //Tetra
00192 Albany::MpasSTKMeshStruct::MpasSTKMeshStruct(const Teuchos::RCP<Teuchos::ParameterList>& params,
00193                                              const Teuchos::RCP<const Epetra_Comm>& comm,
00194                                              const std::vector<int>& indexToTriangleID, int nGlobalTriangles, int numLayers, int Ordering) :
00195   GenericSTKMeshStruct(params,Teuchos::null,3),
00196   out(Teuchos::VerboseObjectBase::getDefaultOStream()),
00197   periodic(false),
00198   NumEles(indexToTriangleID.size()),
00199   hasRestartSol(false),
00200   restartTime(0.)
00201 {
00202   std::vector<int> indexToTetraID(3*indexToTriangleID.size()*numLayers);
00203 
00204   //Int ElemColumnShift = (ordering == ColumnWise) ? 1 : indexToTriangleID.size();
00205   int elemColumnShift = (Ordering == 1) ? 3 : 3*nGlobalTriangles;
00206   int lElemColumnShift = (Ordering == 1) ? 3 : 3*indexToTriangleID.size();
00207   int elemLayerShift = (Ordering == 0) ? 3 : 3*numLayers;
00208 
00209   for(int il=0; il< numLayers; il++)
00210   {
00211     int shift = il*elemColumnShift;
00212     int lShift = il*lElemColumnShift;
00213     for(int j=0; j< indexToTriangleID.size(); j++)
00214     {
00215       for(int iTetra=0; iTetra<3; iTetra++)
00216       {
00217         int lid = lShift + j*elemLayerShift +iTetra;
00218         indexToTetraID[lid] = shift+elemLayerShift * indexToTriangleID[j] +iTetra;
00219       }
00220     }
00221   }
00222 
00223   elem_map = Teuchos::rcp(new Epetra_Map(3*nGlobalTriangles*numLayers, indexToTetraID.size(), &indexToTetraID[0], 0, *comm)); // Distribute the elems equally
00224 
00225   params->validateParameters(*getValidDiscretizationParameters(),0);
00226 
00227 
00228   std::string ebn="Element Block 0";
00229   partVec[0] = & metaData->declare_part(ebn, metaData->element_rank() );
00230   ebNameToIndex[ebn] = 0;
00231 
00232 #ifdef ALBANY_SEACAS
00233   stk::io::put_io_part_attribute(*partVec[0]);
00234 #endif
00235 
00236 
00237   std::vector<std::string> nsNames;
00238   std::string nsn="Lateral";
00239   nsNames.push_back(nsn);
00240   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00241 #ifdef ALBANY_SEACAS
00242     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00243 #endif
00244   nsn="Internal";
00245   nsNames.push_back(nsn);
00246   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00247 #ifdef ALBANY_SEACAS
00248     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00249 #endif
00250   nsn="Bottom";
00251   nsNames.push_back(nsn);
00252   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00253 #ifdef ALBANY_SEACAS
00254   stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00255 #endif
00256 
00257 
00258   std::vector<std::string> ssNames;
00259   std::string ssnLat="lateralside";
00260   std::string ssnBottom="basalside";
00261   std::string ssnTop="upperside";
00262 
00263   ssNames.push_back(ssnLat);
00264   ssNames.push_back(ssnBottom);
00265   ssNames.push_back(ssnTop);
00266   ssPartVec[ssnLat] = & metaData->declare_part(ssnLat, metaData->side_rank() );
00267   ssPartVec[ssnBottom] = & metaData->declare_part(ssnBottom, metaData->side_rank() );
00268   ssPartVec[ssnTop] = & metaData->declare_part(ssnTop, metaData->side_rank() );
00269 #ifdef ALBANY_SEACAS
00270     stk::io::put_io_part_attribute(*ssPartVec[ssnLat]);
00271     stk::io::put_io_part_attribute(*ssPartVec[ssnBottom]);
00272     stk::io::put_io_part_attribute(*ssPartVec[ssnTop]);
00273 #endif
00274 
00275   stk::mesh::fem::set_cell_topology<shards::Tetrahedron<4> >(*partVec[0]);
00276   stk::mesh::fem::set_cell_topology<shards::Triangle<3> >(*ssPartVec[ssnBottom]);
00277   stk::mesh::fem::set_cell_topology<shards::Triangle<3> >(*ssPartVec[ssnTop]);
00278   stk::mesh::fem::set_cell_topology<shards::Triangle<3> >(*ssPartVec[ssnLat]);
00279 
00280   numDim = 3;
00281   int cub = params->get("Cubature Degree",3);
00282   int worksetSizeMax = params->get("Workset Size",50);
00283   int worksetSize = this->computeWorksetSize(worksetSizeMax, elem_map->NumMyElements());
00284 
00285   const CellTopologyData& ctd = *metaData->get_cell_topology(*partVec[0]).getCellTopologyData();
00286 
00287   this->meshSpecs[0] = Teuchos::rcp(new Albany::MeshSpecsStruct(ctd, numDim, cub,
00288                              nsNames, ssNames, worksetSize, partVec[0]->name(),
00289                              ebNameToIndex, this->interleavedOrdering));
00290 
00291 
00292 }
00293 
00294 
00295 Albany::MpasSTKMeshStruct::~MpasSTKMeshStruct()
00296 {
00297 }
00298 
00299 void
00300 Albany::MpasSTKMeshStruct::constructMesh(
00301                                                const Teuchos::RCP<const Epetra_Comm>& comm,
00302                                                const Teuchos::RCP<Teuchos::ParameterList>& params,
00303                                                const unsigned int neq_,
00304                                                const Albany::AbstractFieldContainer::FieldContainerRequirements& req,
00305                                                const Teuchos::RCP<Albany::StateInfoStruct>& sis,
00306                                                const std::vector<int>& indexToVertexID, const std::vector<double>& verticesCoords, const std::vector<bool>& isVertexBoundary, int nGlobalVertices,
00307                                                const std::vector<int>& verticesOnTria,
00308                                                const std::vector<bool>& isBoundaryEdge, const std::vector<int>& trianglesOnEdge, const std::vector<int>& trianglesPositionsOnEdge,
00309                                                const std::vector<int>& verticesOnEdge,
00310                                                const std::vector<int>& indexToEdgeID, int nGlobalEdges,
00311                                                const std::vector<int>& indexToTriangleID,
00312                                                const unsigned int worksetSize,
00313                                                int numLayers, int Ordering)
00314 {
00315   this->SetupFieldData(comm, neq_, req, sis, worksetSize);
00316 
00317     int elemColumnShift = (Ordering == 1) ? 1 : elem_map->NumGlobalElements()/numLayers;
00318     int lElemColumnShift = (Ordering == 1) ? 1 : indexToTriangleID.size();
00319     int elemLayerShift = (Ordering == 0) ? 1 : numLayers;
00320 
00321     int vertexColumnShift = (Ordering == 1) ? 1 : nGlobalVertices;
00322     int lVertexColumnShift = (Ordering == 1) ? 1 : indexToVertexID.size();
00323     int vertexLayerShift = (Ordering == 0) ? 1 : numLayers+1;
00324 
00325     int edgeColumnShift = (Ordering == 1) ? 1 : nGlobalEdges;
00326     int lEdgeColumnShift = (Ordering == 1) ? 1 : indexToEdgeID.size();
00327     int edgeLayerShift = (Ordering == 0) ? 1 : numLayers;
00328 
00329 
00330   metaData->commit();
00331 
00332   bulkData->modification_begin(); // Begin modifying the mesh
00333 
00334   stk::mesh::PartVector nodePartVec;
00335   stk::mesh::PartVector singlePartVec(1);
00336   stk::mesh::PartVector emptyPartVec;
00337   std::cout << "elem_map # elments: " << elem_map->NumMyElements() << std::endl;
00338   unsigned int ebNo = 0; //element block #???
00339 
00340   singlePartVec[0] = nsPartVec["Bottom"];
00341 
00342 
00343   AbstractSTKFieldContainer::IntScalarFieldType* proc_rank_field = fieldContainer->getProcRankField();
00344   AbstractSTKFieldContainer::VectorFieldType* coordinates_field = fieldContainer->getCoordinatesField();
00345   AbstractSTKFieldContainer::ScalarFieldType* surfaceHeight_field = fieldContainer->getSurfaceHeightField();
00346 
00347   for(int i=0; i< (numLayers+1)*indexToVertexID.size(); i++)
00348   {
00349     int ib = (Ordering == 0)*(i%lVertexColumnShift) + (Ordering == 1)*(i/vertexLayerShift);
00350     int il = (Ordering == 0)*(i/lVertexColumnShift) + (Ordering == 1)*(i%vertexLayerShift);
00351 
00352     stk::mesh::Entity* node;
00353     if(il == 0)
00354       node = &bulkData->declare_entity(metaData->node_rank(), il*vertexColumnShift+vertexLayerShift * indexToVertexID[ib]+1, singlePartVec);
00355     else
00356       node = &bulkData->declare_entity(metaData->node_rank(), il*vertexColumnShift+vertexLayerShift * indexToVertexID[ib]+1, nodePartVec);
00357     int numBdEdges(0);
00358     for (int i=0; i<indexToEdgeID.size(); i++)
00359       numBdEdges += isBoundaryEdge[i];
00360 
00361 
00362       double* coord = stk::mesh::field_data(*coordinates_field, *node);
00363     coord[0] = verticesCoords[3*ib];   coord[1] = verticesCoords[3*ib+1]; coord[2] = double(il)/numLayers;
00364 
00365     double* sHeight;
00366      sHeight = stk::mesh::field_data(*surfaceHeight_field, *node);
00367      sHeight[0] = 1.;
00368   }
00369 
00370   for (int i=0; i<elem_map->NumMyElements(); i++) {
00371 
00372    int ib = (Ordering == 0)*(i%lElemColumnShift) + (Ordering == 1)*(i/elemLayerShift);
00373    int il = (Ordering == 0)*(i/lElemColumnShift) + (Ordering == 1)*(i%elemLayerShift);
00374 
00375    int shift = il*vertexColumnShift;
00376 
00377    singlePartVec[0] = partVec[ebNo];
00378      stk::mesh::Entity& elem  = bulkData->declare_entity(metaData->element_rank(), elem_map->GID(i)+1, singlePartVec);
00379 
00380      for(int j=0; j<3; j++)
00381      {
00382        int lowerId = shift+vertexLayerShift * indexToVertexID[verticesOnTria[3*ib+j]]+1;
00383        stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), lowerId);
00384        bulkData->declare_relation(elem, node, j);
00385 
00386        stk::mesh::Entity& node_top = *bulkData->get_entity(metaData->node_rank(), lowerId+vertexColumnShift);
00387        bulkData->declare_relation(elem, node_top, j+3);
00388      }
00389 
00390      int* p_rank = (int*)stk::mesh::field_data(*proc_rank_field, elem);
00391      p_rank[0] = comm->MyPID();
00392   }
00393 
00394 
00395   singlePartVec[0] = ssPartVec["lateralside"];
00396 
00397   //first we store the lateral faces of prisms, which corresponds to edges of the basal mesh
00398 
00399   for (int i=0; i<indexToEdgeID.size()*numLayers; i++) {
00400    int ib = (Ordering == 0)*(i%lEdgeColumnShift) + (Ordering == 1)*(i/edgeLayerShift);
00401    if(isBoundaryEdge[ib])
00402    {
00403      int il = (Ordering == 0)*(i/lEdgeColumnShift) + (Ordering == 1)*(i%edgeLayerShift);
00404      int basalEdgeId = indexToEdgeID[ib]*edgeLayerShift;
00405      int basalElemId = indexToTriangleID[trianglesOnEdge[2*ib]]*elemLayerShift;
00406      int basalVertexId[2] = {indexToVertexID[verticesOnEdge[2*ib]]*vertexLayerShift, indexToVertexID[verticesOnEdge[2*ib+1]]*vertexLayerShift};
00407      stk::mesh::Entity& side = bulkData->declare_entity(metaData->side_rank(), edgeColumnShift*il+basalEdgeId+1, singlePartVec);
00408      stk::mesh::Entity& elem  = *bulkData->get_entity(metaData->element_rank(),  basalElemId+elemColumnShift*il+1);
00409      bulkData->declare_relation(elem, side,  trianglesPositionsOnEdge[2*ib] );
00410      for(int j=0; j<2; j++)
00411      {
00412        stk::mesh::Entity& nodeBottom = *bulkData->get_entity(metaData->node_rank(), basalVertexId[j]+vertexColumnShift*il+1);
00413        bulkData->declare_relation(side, nodeBottom, j);
00414        stk::mesh::Entity& nodeTop = *bulkData->get_entity(metaData->node_rank(), basalVertexId[j]+vertexColumnShift*(il+1)+1);
00415        bulkData->declare_relation(side, nodeTop, j+2);
00416      }
00417    }
00418   }
00419 
00420   //then we store the lower and upper faces of prisms, which corresponds to triangles of the basal mesh
00421 
00422   edgeLayerShift = (Ordering == 0) ? 1 : numLayers+1;
00423   edgeColumnShift = elemColumnShift;
00424 
00425   singlePartVec[0] = ssPartVec["basalside"];
00426 
00427   int edgeOffset = nGlobalEdges*numLayers;
00428   for (int i=0; i<indexToTriangleID.size(); i++)
00429   {
00430     stk::mesh::Entity& side = bulkData->declare_entity(metaData->side_rank(), indexToTriangleID[i]*edgeLayerShift+edgeOffset+1, singlePartVec);
00431     stk::mesh::Entity& elem  = *bulkData->get_entity(metaData->element_rank(),  indexToTriangleID[i]*elemLayerShift+1);
00432     bulkData->declare_relation(elem, side,  3);
00433     for(int j=0; j<3; j++)
00434     {
00435      stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), vertexLayerShift*indexToVertexID[verticesOnTria[3*i+j]]+1);
00436      bulkData->declare_relation(side, node, j);
00437     }
00438   }
00439 
00440   singlePartVec[0] = ssPartVec["upperside"];
00441 
00442   for (int i=0; i<indexToTriangleID.size(); i++)
00443   {
00444       stk::mesh::Entity& side = bulkData->declare_entity(metaData->side_rank(), indexToTriangleID[i]*edgeLayerShift+numLayers*edgeColumnShift+edgeOffset+1, singlePartVec);
00445       stk::mesh::Entity& elem  = *bulkData->get_entity(metaData->element_rank(),  indexToTriangleID[i]*elemLayerShift+(numLayers-1)*elemColumnShift+1);
00446       bulkData->declare_relation(elem, side,  4);
00447       for(int j=0; j<3; j++)
00448       {
00449        stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), vertexLayerShift*indexToVertexID[verticesOnTria[3*i+j]]+numLayers*vertexColumnShift+1);
00450        bulkData->declare_relation(side, node, j);
00451       }
00452   }
00453 
00454   bulkData->modification_end();
00455 }
00456 
00457 
00458 
00459 void
00460 Albany::MpasSTKMeshStruct::constructMesh(
00461                                                const Teuchos::RCP<const Epetra_Comm>& comm,
00462                                                const Teuchos::RCP<Teuchos::ParameterList>& params,
00463                                                const unsigned int neq_,
00464                                                const Albany::AbstractFieldContainer::FieldContainerRequirements& req,
00465                                                const Teuchos::RCP<Albany::StateInfoStruct>& sis,
00466                                                const std::vector<int>& indexToVertexID, const std::vector<int>& indexToMpasVertexID, const std::vector<double>& verticesCoords, const std::vector<bool>& isVertexBoundary, int nGlobalVertices,
00467                                                const std::vector<int>& verticesOnTria,
00468                                                const std::vector<bool>& isBoundaryEdge, const std::vector<int>& trianglesOnEdge, const std::vector<int>& trianglesPositionsOnEdge,
00469                                                const std::vector<int>& verticesOnEdge,
00470                                                const std::vector<int>& indexToEdgeID, int nGlobalEdges,
00471                                                const std::vector<int>& indexToTriangleID,
00472                                                const unsigned int worksetSize,
00473                                                int numLayers, int Ordering)
00474 {
00475   this->SetupFieldData(comm, neq_, req, sis, worksetSize);
00476 
00477     int elemColumnShift = (Ordering == 1) ? 3 : elem_map->NumGlobalElements()/numLayers;
00478     int lElemColumnShift = (Ordering == 1) ? 3 : 3*indexToTriangleID.size();
00479     int elemLayerShift = (Ordering == 0) ? 3 : 3*numLayers;
00480 
00481     int vertexColumnShift = (Ordering == 1) ? 1 : nGlobalVertices;
00482     int lVertexColumnShift = (Ordering == 1) ? 1 : indexToVertexID.size();
00483     int vertexLayerShift = (Ordering == 0) ? 1 : numLayers+1;
00484 
00485     int edgeColumnShift = (Ordering == 1) ? 2 : 2*nGlobalEdges;
00486     int lEdgeColumnShift = (Ordering == 1) ? 1 : indexToEdgeID.size();
00487     int edgeLayerShift = (Ordering == 0) ? 1 : numLayers;
00488 
00489 
00490   metaData->commit();
00491 
00492   bulkData->modification_begin(); // Begin modifying the mesh
00493 
00494   stk::mesh::PartVector nodePartVec;
00495   stk::mesh::PartVector singlePartVec(1);
00496   stk::mesh::PartVector emptyPartVec;
00497   std::cout << "elem_map # elments: " << elem_map->NumMyElements() << std::endl;
00498   unsigned int ebNo = 0; //element block #???
00499 
00500   singlePartVec[0] = nsPartVec["Bottom"];
00501 
00502   AbstractSTKFieldContainer::IntScalarFieldType* proc_rank_field = fieldContainer->getProcRankField();
00503   AbstractSTKFieldContainer::VectorFieldType* coordinates_field = fieldContainer->getCoordinatesField();
00504   AbstractSTKFieldContainer::ScalarFieldType* surfaceHeight_field = fieldContainer->getSurfaceHeightField();
00505 
00506 
00507   for(int i=0; i< (numLayers+1)*indexToVertexID.size(); i++)
00508   {
00509     int ib = (Ordering == 0)*(i%lVertexColumnShift) + (Ordering == 1)*(i/vertexLayerShift);
00510     int il = (Ordering == 0)*(i/lVertexColumnShift) + (Ordering == 1)*(i%vertexLayerShift);
00511 
00512     stk::mesh::Entity* node;
00513     if(il == 0)
00514       node = &bulkData->declare_entity(metaData->node_rank(), il*vertexColumnShift+vertexLayerShift * indexToVertexID[ib]+1, singlePartVec);
00515     else
00516       node = &bulkData->declare_entity(metaData->node_rank(), il*vertexColumnShift+vertexLayerShift * indexToVertexID[ib]+1, nodePartVec);
00517 
00518       double* coord = stk::mesh::field_data(*coordinates_field, *node);
00519     coord[0] = verticesCoords[3*ib];   coord[1] = verticesCoords[3*ib+1]; coord[2] = double(il)/numLayers;
00520 
00521     double* sHeight;
00522      sHeight = stk::mesh::field_data(*surfaceHeight_field, *node);
00523      sHeight[0] = 1.;
00524   }
00525 
00526   int tetrasLocalIdsOnPrism[3][4];
00527 
00528   for (int i=0; i<elem_map->NumMyElements()/3; i++) {
00529 
00530    int ib = (Ordering == 0)*(i%(lElemColumnShift/3)) + (Ordering == 1)*(i/(elemLayerShift/3));
00531    int il = (Ordering == 0)*(i/(lElemColumnShift/3)) + (Ordering == 1)*(i%(elemLayerShift/3));
00532 
00533    int shift = il*vertexColumnShift;
00534 
00535    singlePartVec[0] = partVec[ebNo];
00536 
00537 
00538      //TODO: this could be done only in the first layer and then copied into the other layers
00539      int prismMpasIds[3], prismGlobalIds[6];
00540      for (int j = 0; j < 3; j++)
00541    {
00542        int mpasLowerId = vertexLayerShift * indexToMpasVertexID[verticesOnTria[3*ib+j]];
00543        int lowerId = shift+vertexLayerShift * indexToVertexID[verticesOnTria[3*ib+j]];
00544        prismMpasIds[j] = mpasLowerId;
00545      prismGlobalIds[j] = lowerId;
00546      prismGlobalIds[j + 3] = lowerId+vertexColumnShift;
00547    }
00548 
00549      tetrasFromPrismStructured (prismMpasIds, prismGlobalIds, tetrasLocalIdsOnPrism);
00550 
00551 
00552      for(int iTetra = 0; iTetra<3; iTetra++)
00553      {
00554        stk::mesh::Entity& elem  = bulkData->declare_entity(metaData->element_rank(), elem_map->GID(3*i+iTetra)+1, singlePartVec);
00555      for(int j=0; j<4; j++)
00556      {
00557        stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), tetrasLocalIdsOnPrism[iTetra][j]+1);
00558        bulkData->declare_relation(elem, node, j);
00559      }
00560      int* p_rank = (int*)stk::mesh::field_data(*proc_rank_field, elem);
00561      p_rank[0] = comm->MyPID();
00562      }
00563 
00564 
00565   }
00566 
00567 
00568   singlePartVec[0] = ssPartVec["lateralside"];
00569 
00570   //first we store the lateral faces of prisms, which corresponds to edges of the basal mesh
00571   int tetraSidePoints[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
00572   std::vector<int> tetraPos(2), facePos(2);
00573 
00574   std::vector<std::vector<std::vector<int> > > prismStruct(3, std::vector<std::vector<int> >(4, std::vector<int>(3)));
00575   for (int i=0; i<indexToEdgeID.size()*numLayers; i++) {
00576    int ib = (Ordering == 0)*(i%lEdgeColumnShift) + (Ordering == 1)*(i/edgeLayerShift);
00577    if(isBoundaryEdge[ib])
00578    {
00579      int il = (Ordering == 0)*(i/lEdgeColumnShift) + (Ordering == 1)*(i%edgeLayerShift);
00580      int lBasalElemId = trianglesOnEdge[2*ib];
00581      int basalElemId = indexToTriangleID[lBasalElemId];
00582 
00583      //TODO: this could be done only in the first layer and then copied into the other layers
00584      int prismMpasIds[3], prismGlobalIds[6];
00585      int shift = il*vertexColumnShift;
00586      for (int j = 0; j < 3; j++)
00587      {
00588        int mpasLowerId = vertexLayerShift * indexToMpasVertexID[verticesOnTria[3*lBasalElemId+j]];
00589        int lowerId = shift+vertexLayerShift * indexToVertexID[verticesOnTria[3*lBasalElemId+j]];
00590        prismMpasIds[j] = mpasLowerId;
00591        prismGlobalIds[j] = lowerId;
00592        prismGlobalIds[j + 3] = lowerId+vertexColumnShift;
00593      }
00594 
00595       tetrasFromPrismStructured (prismMpasIds, prismGlobalIds, tetrasLocalIdsOnPrism);
00596 
00597 
00598     for(int iTetra = 0; iTetra<3; iTetra++)
00599       {
00600        std::vector<std::vector<int> >& tetraStruct =prismStruct[iTetra];
00601        stk::mesh::EntityId tetraPoints[4];
00602        for(int j=0; j<4; j++)
00603        {
00604          tetraPoints[j] = bulkData->get_entity(metaData->node_rank(), tetrasLocalIdsOnPrism[iTetra][j]+1)->identifier();
00605         // std::cout<< tetraPoints[j] << ", ";
00606        }
00607        for(int iFace=0; iFace<4; iFace++)
00608        {
00609          std::vector<int>&  face = tetraStruct[iFace];
00610          for(int j=0; j<3; j++)
00611            face[j] = tetraPoints[tetraSidePoints[iFace][j]];
00612        }
00613       }
00614 
00615 
00616 
00617      int basalVertexId[2] = {indexToVertexID[verticesOnEdge[2*ib]]*vertexLayerShift, indexToVertexID[verticesOnEdge[2*ib+1]]*vertexLayerShift};
00618      std::vector<int> bdPrismFaceIds(4);
00619 
00620      bdPrismFaceIds[0] = indexToVertexID[verticesOnEdge[2*ib]]*vertexLayerShift+vertexColumnShift*il+1;
00621      bdPrismFaceIds[1] = indexToVertexID[verticesOnEdge[2*ib+1]]*vertexLayerShift+vertexColumnShift*il+1;
00622      bdPrismFaceIds[2] = bdPrismFaceIds[0]+vertexColumnShift;
00623      bdPrismFaceIds[3] = bdPrismFaceIds[1]+vertexColumnShift;
00624 
00625      //std::cout<< "bdPrismFaceIds: (" << bdPrismFaceIds[0] << ", " << bdPrismFaceIds[1] << ", " << bdPrismFaceIds[2] << ", " << bdPrismFaceIds[3] << ")"<<std::endl;
00626 
00627 
00628 
00629      setBdFacesOnPrism (prismStruct, bdPrismFaceIds, tetraPos, facePos);
00630 
00631      int basalEdgeId = indexToEdgeID[ib]*2*edgeLayerShift;
00632      for(int k=0; k< tetraPos.size(); k++)
00633      {
00634        int iTetra = tetraPos[k];
00635        int iFace = facePos[k];
00636        stk::mesh::Entity& elem = *bulkData->get_entity(metaData->element_rank(), il*elemColumnShift+elemLayerShift * basalElemId +iTetra+1);
00637        std::vector<int>& faceIds = prismStruct[iTetra][iFace];
00638        stk::mesh::Entity& side = bulkData->declare_entity(metaData->side_rank(), edgeColumnShift*il+basalEdgeId+k+1, singlePartVec);
00639        bulkData->declare_relation(elem, side,  iFace );
00640        for(int j=0; j<3; j++)
00641        {
00642          stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), faceIds[j]);
00643          bulkData->declare_relation(side, node, j);
00644        }
00645      }
00646    }
00647   }
00648 
00649   //then we store the lower and upper faces of prisms, which corresponds to triangles of the basal mesh
00650 
00651   edgeLayerShift = (Ordering == 0) ? 1 : numLayers+1;
00652   edgeColumnShift = 2*(elemColumnShift/3);
00653 
00654   singlePartVec[0] = ssPartVec["basalside"];
00655 
00656   int edgeOffset = 2*nGlobalEdges*numLayers;
00657   for (int i=0; i<indexToTriangleID.size(); i++)
00658   {
00659     stk::mesh::Entity& side = bulkData->declare_entity(metaData->side_rank(), indexToTriangleID[i]*edgeLayerShift+edgeOffset+1, singlePartVec);
00660     stk::mesh::Entity& elem  = *bulkData->get_entity(metaData->element_rank(),  indexToTriangleID[i]*elemLayerShift+1);
00661     bulkData->declare_relation(elem, side,  3);
00662     for(int j=0; j<3; j++)
00663     {
00664      stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), vertexLayerShift*indexToVertexID[verticesOnTria[3*i+j]]+1);
00665      bulkData->declare_relation(side, node, j);
00666     }
00667   }
00668 
00669   singlePartVec[0] = ssPartVec["upperside"];
00670 
00671   for (int i=0; i<indexToTriangleID.size(); i++)
00672   {
00673   stk::mesh::Entity& side = bulkData->declare_entity(metaData->side_rank(), indexToTriangleID[i]*edgeLayerShift+numLayers*edgeColumnShift+edgeOffset+1, singlePartVec);
00674   stk::mesh::Entity& elem  = *bulkData->get_entity(metaData->element_rank(),  indexToTriangleID[i]*elemLayerShift+(numLayers-1)*elemColumnShift+1+2);
00675   bulkData->declare_relation(elem, side,  1);
00676   for(int j=0; j<3; j++)
00677   {
00678     stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), vertexLayerShift*indexToVertexID[verticesOnTria[3*i+j]]+numLayers*vertexColumnShift+1);
00679     bulkData->declare_relation(side, node, j);
00680   }
00681   }
00682 
00683   bulkData->modification_end();
00684 }
00685 
00686 void
00687 Albany::MpasSTKMeshStruct::constructMesh(
00688                                                const Teuchos::RCP<const Epetra_Comm>& comm,
00689                                                const Teuchos::RCP<Teuchos::ParameterList>& params,
00690                                                const unsigned int neq_,
00691                                                const Albany::AbstractFieldContainer::FieldContainerRequirements& req,
00692                                                const Teuchos::RCP<Albany::StateInfoStruct>& sis,
00693                                                const std::vector<int>& indexToVertexID, const std::vector<double>& verticesCoords, const std::vector<bool>& isVertexBoundary, int nGlobalVertices,
00694                                                const std::vector<int>& verticesOnTria,
00695                                                const std::vector<bool>& isBoundaryEdge, const std::vector<int>& trianglesOnEdge, const std::vector<int>& trianglesPositionsOnEdge,
00696                                                const std::vector<int>& verticesOnEdge,
00697                                                const std::vector<int>& indexToEdgeID, int nGlobalEdges,
00698                                                const unsigned int worksetSize)
00699 {
00700   this->SetupFieldData(comm, neq_, req, sis, worksetSize);
00701 
00702   metaData->commit();
00703 
00704   bulkData->modification_begin(); // Begin modifying the mesh
00705 
00706   stk::mesh::PartVector nodePartVec;
00707   stk::mesh::PartVector singlePartVec(1);
00708   stk::mesh::PartVector emptyPartVec;
00709   std::cout << "elem_map # elments: " << elem_map->NumMyElements() << std::endl;
00710   unsigned int ebNo = 0; //element block #??? 
00711   int sideID = 0;
00712 
00713   AbstractSTKFieldContainer::IntScalarFieldType* proc_rank_field = fieldContainer->getProcRankField();
00714   AbstractSTKFieldContainer::VectorFieldType* coordinates_field = fieldContainer->getCoordinatesField();
00715   AbstractSTKFieldContainer::ScalarFieldType* surfaceHeight_field = fieldContainer->getSurfaceHeightField();
00716 
00717   for (int i=0; i<indexToVertexID.size(); i++)
00718   {
00719     stk::mesh::Entity& node = bulkData->declare_entity(metaData->node_rank(), indexToVertexID[i]+1, nodePartVec);
00720 
00721     double* coord;
00722     coord = stk::mesh::field_data(*coordinates_field, node);
00723     coord[0] = verticesCoords[3*i];   coord[1] = verticesCoords[3*i+1]; coord[2] = verticesCoords[3*i+2];
00724   }
00725 
00726   for (int i=0; i<elem_map->NumMyElements(); i++)
00727   {
00728 
00729      singlePartVec[0] = partVec[ebNo];
00730      stk::mesh::Entity& elem  = bulkData->declare_entity(metaData->element_rank(), elem_map->GID(i)+1, singlePartVec);
00731 
00732      for(int j=0; j<3; j++)
00733      {
00734        stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), indexToVertexID[verticesOnTria[3*i+j]]+1);
00735        bulkData->declare_relation(elem, node, j);
00736      }
00737     
00738      int* p_rank = (int*)stk::mesh::field_data(*proc_rank_field, elem);
00739      p_rank[0] = comm->MyPID();
00740   }
00741 
00742   for (int i=0; i<indexToEdgeID.size(); i++) {
00743 
00744    if(isBoundaryEdge[i])
00745    {
00746 
00747      singlePartVec[0] = ssPartVec["lateralside"];
00748      stk::mesh::Entity& side = bulkData->declare_entity(metaData->side_rank(), indexToEdgeID[i]+1, singlePartVec);
00749      stk::mesh::Entity& elem  = *bulkData->get_entity(metaData->element_rank(),  elem_map->GID(trianglesOnEdge[2*i])+1);
00750      bulkData->declare_relation(elem, side,  trianglesPositionsOnEdge[2*i] /*local side id*/);
00751      for(int j=0; j<2; j++)
00752      {
00753        stk::mesh::Entity& node = *bulkData->get_entity(metaData->node_rank(), indexToVertexID[verticesOnEdge[2*i+j]]+1);
00754        bulkData->declare_relation(side, node, j);
00755      }
00756    }
00757   }
00758 
00759   bulkData->modification_end();
00760 }
00761 
00762 Teuchos::RCP<const Teuchos::ParameterList>
00763 Albany::MpasSTKMeshStruct::getValidDiscretizationParameters() const
00764 {
00765   Teuchos::RCP<Teuchos::ParameterList> validPL =
00766     this->getValidGenericSTKParameters("Valid ASCII_DiscParams");
00767 
00768   return validPL;
00769 }