Albany: a Trilinos-based PDE code: Albany_AsciiSTKMeshStruct.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_AsciiSTKMeshStruct.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 
00020 #ifdef ALBANY_SEACAS
00021 #include <stk_io/IossBridge.hpp>
00022 #endif
00023 
00024 
00025 //#include <stk_mesh/fem/FEMHelpers.hpp>
00026 #include <boost/algorithm/string/predicate.hpp>
00027 
00028 #include "Albany_Utils.hpp"
00029 
00030 
00031 //uncomment the following line if you want debug output to be printed to screen
00032 //#define OUTPUT_TO_SCREEN
00033 
00034 
00035 //Constructor for meshes read from ASCII file 
00036 Albany::AsciiSTKMeshStruct::AsciiSTKMeshStruct(
00037                                              const Teuchos::RCP<Teuchos::ParameterList>& params, 
00038                                              const Teuchos::RCP<const Epetra_Comm>& comm) :
00039   GenericSTKMeshStruct(params,Teuchos::null,3),
00040   out(Teuchos::VerboseObjectBase::getDefaultOStream()),
00041   periodic(false)
00042 {
00043    int numProc = comm->NumProc(); //total number of processors
00044    contigIDs = params->get("Contiguous IDs", true); 
00045    std::cout << "Number of processors: " << numProc << std::endl; 
00046    //names of files giving the mesh
00047    char meshfilename[100]; 
00048    char shfilename[100];
00049    char confilename[100];
00050    char bffilename[100];
00051    char geIDsfilename[100];
00052    char gnIDsfilename[100];
00053    char bfIDsfilename[100];
00054    char flwafilename[100]; //flow factor file
00055    char tempfilename[100]; //temperature file
00056    char betafilename[100]; //basal friction coefficient file
00057    if ((numProc == 1) & (contigIDs == true)) { //serial run with contiguous global IDs
00058 #ifdef OUTPUT_TO_SCREEN
00059      std::cout << "Ascii mesh has contiguous IDs; no bfIDs, geIDs, gnIDs files required." << std::endl;
00060 #endif
00061      sprintf(meshfilename, "%s", "xyz");
00062      sprintf(shfilename, "%s", "sh");
00063      sprintf(confilename, "%s", "eles");
00064      sprintf(bffilename, "%s", "bf");
00065      sprintf(flwafilename, "%s", "flwa");
00066      sprintf(tempfilename, "%s", "temp");
00067      sprintf(betafilename, "%s", "beta");
00068    }
00069    else { //parallel run or serial run with non-contiguous global IDs - proc # is appended to file name to indicate what processor the mesh piece is on 
00070 #ifdef OUTPUT_TO_SCREEN
00071      if ((numProc == 1) & (contigIDs == false))
00072         std::cout << "1 processor run with non-contiguous IDs; bfIDs0, geIDs0, gnIDs0 files required." << std::endl;
00073 #endif
00074      int suffix = comm->MyPID(); //current processor number 
00075      sprintf(meshfilename, "%s%i", "xyz", suffix);
00076      sprintf(shfilename, "%s%i", "sh", suffix);
00077      sprintf(confilename, "%s%i", "eles", suffix);
00078      sprintf(bffilename, "%s%i", "bf", suffix);
00079      sprintf(geIDsfilename, "%s%i", "geIDs", suffix);
00080      sprintf(gnIDsfilename, "%s%i", "gnIDs", suffix);
00081      sprintf(bfIDsfilename, "%s%i", "bfIDs", suffix);
00082      sprintf(flwafilename, "%s%i", "flwa", suffix);
00083      sprintf(tempfilename, "%s%i", "temp", suffix);
00084      sprintf(betafilename, "%s%i", "beta", suffix);
00085    }
00086 
00087     //read in coordinates of mesh -- right now hard coded for 3D
00088     //assumes mesh file is called "xyz" and its first row is the number of nodes  
00089     FILE *meshfile = fopen(meshfilename,"r");
00090     if (meshfile == NULL) { //check if coordinates file exists
00091       *out << "Error in AsciiSTKMeshStruct: coordinates file " << meshfilename <<" not found!"<< std::endl;
00092       TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00093           std::endl << "Error in AsciiSTKMeshStruct: coordinates file " << meshfilename << " not found!"<< std::endl);
00094     }
00095     double temp; 
00096     fseek(meshfile, 0, SEEK_SET); 
00097     fscanf(meshfile, "%lf", &temp); 
00098     NumNodes = int(temp); 
00099 #ifdef OUTPUT_TO_SCREEN
00100     *out << "numNodes: " << NumNodes << std::endl;  
00101 #endif
00102     xyz = new double[NumNodes][3]; 
00103     char buffer[100];
00104     fgets(buffer, 100, meshfile); 
00105     for (int i=0; i<NumNodes; i++){
00106       fgets(buffer, 100, meshfile); 
00107       sscanf(buffer, "%lf %lf %lf", &xyz[i][0], &xyz[i][1], &xyz[i][2]); 
00108       //*out << "i: " << i << ", x: " << xyz[i][0] << ", y: " << xyz[i][1] << ", z: " << xyz[i][2] << std::endl; 
00109      }
00110     //read in surface height data from mesh 
00111     //assumes surface height file is called "sh" and its first row is the number of nodes  
00112     FILE *shfile = fopen(shfilename,"r");
00113     have_sh = false;
00114     if (shfile != NULL) have_sh = true;
00115     if (have_sh) {
00116       fseek(shfile, 0, SEEK_SET); 
00117       fscanf(shfile, "%lf", &temp); 
00118       int NumNodesSh = int(temp);
00119 #ifdef OUTPUT_TO_SCREEN
00120       *out << "NumNodesSh: " << NumNodesSh<< std::endl;
00121 #endif 
00122       if (NumNodesSh != NumNodes) { 
00123            *out << "Error in AsciiSTKMeshStruct: sh file must have same number nodes as xyz file!  numNodes in xyz = " << NumNodes <<", numNodes in sh = "<< NumNodesSh  << std::endl;
00124           TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00125             std::endl << "Error in AsciiSTKMeshStruct: sh file must have same number nodes as xyz file!  numNodes in xyz = " << NumNodes << ", numNodes in sh = "<< NumNodesSh << std::endl);
00126       }
00127       sh = new double[NumNodes]; 
00128       fgets(buffer, 100, shfile); 
00129       for (int i=0; i<NumNodes; i++){
00130         fgets(buffer, 100, shfile); 
00131         sscanf(buffer, "%lf", &sh[i]); 
00132         //*out << "i: " << i << ", sh: " << sh[i] << std::endl; 
00133        }
00134      }
00135      //read in connectivity file -- right now hard coded for 3D hexes
00136      //assumes mesh file is called "eles" and its first row is the number of elements  
00137      FILE *confile = fopen(confilename,"r"); 
00138      if (confile == NULL) { //check if element connectivity file exists
00139       *out << "Error in AsciiSTKMeshStruct: element connectivity file " << confilename <<" not found!"<< std::endl;
00140       TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00141           std::endl << "Error in AsciiSTKMeshStruct: element connectivity file " << confilename << " not found!"<< std::endl);
00142      }
00143      fseek(confile, 0, SEEK_SET); 
00144      fscanf(confile, "%lf", &temp); 
00145      NumEles = int(temp); 
00146 #ifdef OUTPUT_TO_SCREEN
00147      *out << "numEles: " << NumEles << std::endl;
00148 #endif 
00149      eles = new int[NumEles][8]; 
00150      fgets(buffer, 100, confile); 
00151      for (int i=0; i<NumEles; i++){
00152         fgets(buffer, 100, confile); 
00153         sscanf(buffer, "%i %i %i %i %i %i %i %i", &eles[i][0], &eles[i][1], &eles[i][2], &eles[i][3], &eles[i][4], &eles[i][5], &eles[i][6], &eles[i][7]);
00154         //*out << "elt # " << i << ": " << eles[i][0] << " " << eles[i][1] << " " << eles[i][2] << " " << eles[i][3] << " " << eles[i][4] << " "
00155         //                  << eles[i][5] << " " << eles[i][6] << " " << eles[i][7] << std::endl; 
00156      }
00157     //read in basal face connectivity file from ascii file
00158     //assumes basal face connectivity file is called "bf" and its first row is the number of faces on basal boundary
00159     FILE *bffile = fopen(bffilename,"r");
00160     have_bf = false;
00161     if (bffile != NULL) have_bf = true;
00162     if (have_bf) {
00163       fseek(bffile, 0, SEEK_SET); 
00164       fscanf(bffile, "%lf", &temp); 
00165       NumBasalFaces = int(temp); 
00166 #ifdef OUTPUT_TO_SCREEN
00167       *out << "numBasalFaces: " << NumBasalFaces << std::endl;  
00168 #endif
00169       bf = new int[NumBasalFaces][5]; //1st column of bf: element # that face belongs to, 2rd-5th columns of bf: connectivity (hard-coded for quad faces) 
00170       fgets(buffer, 100, bffile); 
00171       for (int i=0; i<NumBasalFaces; i++){
00172         fgets(buffer, 100, bffile); 
00173         sscanf(buffer, "%i %i %i %i %i", &bf[i][0], &bf[i][1], &bf[i][2], &bf[i][3], &bf[i][4]); 
00174         //*out << "face #:" << bf[i][0] << ", face conn:" << bf[i][1] << " " << bf[i][2] << " " << bf[i][3] << " " << bf[i][4] << std::endl; 
00175        }
00176      }
00177      //Create array w/ global element IDs 
00178      globalElesID = new int[NumEles];
00179      if ((numProc == 1) & (contigIDs == true)) { //serial run with contiguous global IDs: element IDs are just 0->NumEles-1
00180        for (int i=0; i<NumEles; i++) {
00181           globalElesID[i] = i; 
00182           //*out << "local element ID #:" << i << ", global element ID #:" << globalElesID[i] << std::endl;
00183        }
00184      }
00185      else {//parallel run: read global element IDs from file.  
00186            //This file should have a header like the other files, and length NumEles.
00187        FILE *geIDsfile = fopen(geIDsfilename,"r");
00188        if (geIDsfile == NULL) { //check if global element IDs file exists
00189          *out << "Error in AsciiSTKMeshStruct: global element IDs file " << geIDsfilename <<" not found!"<< std::endl;
00190          TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00191             std::endl << "Error in AsciiSTKMeshStruct: global element IDs file " << geIDsfilename << " not found!"<< std::endl);
00192        }
00193        fseek(geIDsfile, 0, SEEK_SET);
00194        fgets(buffer, 100, geIDsfile);
00195        for (int i=0; i<NumEles; i++){
00196          fgets(buffer, 100, geIDsfile);
00197          sscanf(buffer, "%i ", &globalElesID[i]);
00198          globalElesID[i] = globalElesID[i]-1; //subtract 1 b/c global element IDs file assumed to be 1-based not 0-based
00199          //*out << "local element ID #:" << i << ", global element ID #:" << globalElesID[i] << std::endl;
00200        }
00201      }
00202      //Create array w/ global node IDs 
00203      globalNodesID = new int[NumNodes];
00204      if ((numProc == 1) & (contigIDs == true)) { //serial run with contiguous global IDs: element IDs are just 0->NumEles-1
00205        for (int i=0; i<NumNodes; i++) { 
00206           globalNodesID[i] = i; 
00207           //*out << "local node ID #:" << i << ", global node ID #:" << globalNodesID[i] << std::endl;
00208        }
00209      }
00210      else {//parallel run: read global node IDs from file.  
00211            //This file should have a header like the other files, and length NumNodes
00212        FILE *gnIDsfile = fopen(gnIDsfilename,"r");
00213        if (gnIDsfile == NULL) { //check if global node IDs file exists
00214          *out << "Error in AsciiSTKMeshStruct: global node IDs file " << gnIDsfilename <<" not found!"<< std::endl;
00215          TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00216             std::endl << "Error in AsciiSTKMeshStruct: global node IDs file " << gnIDsfilename << " not found!"<< std::endl);
00217        }
00218        fseek(gnIDsfile, 0, SEEK_SET);
00219        fgets(buffer, 100, gnIDsfile);
00220        for (int i=0; i<NumNodes; i++){
00221          fgets(buffer, 100, gnIDsfile);
00222          sscanf(buffer, "%i ", &globalNodesID[i]);
00223          globalNodesID[i] = globalNodesID[i]-1; //subtract 1 b/c global node IDs file assumed to be 1-based not 0-based 
00224          //*out << "local node ID #:" << i << ", global node ID #:" << globalNodesID[i] << std::endl;
00225        }
00226      }
00227      basalFacesID = new int[NumBasalFaces];
00228      if ((numProc == 1) & (contigIDs == true)) { //serial run with contiguous global IDs: element IDs are just 0->NumEles-1
00229        for (int i=0; i<NumBasalFaces; i++) { 
00230           basalFacesID[i] = i; 
00231           //*out << "local face ID #:" << i << ", global face ID #:" << basalFacesID[i] << std::endl;
00232        }
00233      }
00234      else {//parallel run: read basal face IDs from file.  
00235            //This file should have a header like the other files, and length NumBasalFaces
00236        FILE *bfIDsfile = fopen(bfIDsfilename,"r");
00237        fseek(bfIDsfile, 0, SEEK_SET);
00238        fgets(buffer, 100, bfIDsfile);
00239        for (int i=0; i<NumBasalFaces; i++){
00240          fgets(buffer, 100, bfIDsfile);
00241          sscanf(buffer, "%i ", &basalFacesID[i]);
00242          basalFacesID[i] = basalFacesID[i]-1; //subtract 1 b/c basal face IDs file assumed to be 1-based not 0-based
00243          //*out << "local face ID #:" << i << ", global face ID #:" << basalFacesID[i] << std::endl;
00244        }
00245      }
00246     //read in flow factor (flwa) data from mesh 
00247     //assumes flow factor file is called "flwa" and its first row is the number of elements in the mesh
00248     FILE *flwafile = fopen(flwafilename,"r");
00249     have_flwa = false;
00250     if (flwafile != NULL) have_flwa = true;
00251     if (have_flwa) {
00252       fseek(flwafile, 0, SEEK_SET); 
00253       fscanf(flwafile, "%lf", &temp); 
00254       flwa = new double[NumEles]; 
00255       fgets(buffer, 100, flwafile); 
00256       for (int i=0; i<NumEles; i++){
00257         fgets(buffer, 100, flwafile); 
00258         sscanf(buffer, "%lf", &flwa[i]); 
00259         //*out << "i: " << i << ", flwa: " << flwa[i] << std::endl; 
00260        }
00261      }
00262     //read in temperature data from mesh 
00263     //assumes temperature file is called "temp" and its first row is the number of elements in the mesh
00264     FILE *tempfile = fopen(tempfilename,"r");
00265     have_temp = false;
00266     if (tempfile != NULL) have_temp = true;
00267     if (have_temp) {
00268       fseek(tempfile, 0, SEEK_SET); 
00269       fscanf(tempfile, "%lf", &temp); 
00270       temper = new double[NumEles]; 
00271       fgets(buffer, 100, tempfile); 
00272       for (int i=0; i<NumEles; i++){
00273         fgets(buffer, 100, tempfile); 
00274         sscanf(buffer, "%lf", &temper[i]); 
00275         //*out << "i: " << i << ", temp: " << temper[i] << std::endl; 
00276        }
00277      }
00278     //read in basal friction (beta) data from mesh 
00279     //assumes basal friction file is called "beta" and its first row is the number of nodes  
00280     FILE *betafile = fopen(betafilename,"r");
00281     have_beta = false;
00282     if (betafile != NULL) have_beta = true;
00283     if (have_beta) {
00284       fseek(betafile, 0, SEEK_SET); 
00285       fscanf(betafile, "%lf", &temp); 
00286       beta = new double[NumNodes]; 
00287       fgets(buffer, 100, betafile); 
00288       for (int i=0; i<NumNodes; i++){
00289         fgets(buffer, 100, betafile); 
00290         sscanf(buffer, "%lf", &beta[i]); 
00291         //*out << "i: " << i << ", beta: " << beta[i] << std::endl; 
00292        }
00293      }
00294  
00295   elem_map = Teuchos::rcp(new Epetra_Map(-1, NumEles, globalElesID, 0, *comm)); //Distribute the elements according to the global element IDs
00296   node_map = Teuchos::rcp(new Epetra_Map(-1, NumNodes, globalNodesID, 0, *comm)); //Distribute the nodes according to the global node IDs 
00297   basal_face_map = Teuchos::rcp(new Epetra_Map(-1, NumBasalFaces, basalFacesID, 0, *comm)); //Distribute the elements according to the basal face IDs
00298 
00299   
00300   params->validateParameters(*getValidDiscretizationParameters(),0);
00301 
00302 
00303   std::string ebn="Element Block 0";
00304   partVec[0] = & metaData->declare_part(ebn, metaData->element_rank() );
00305   ebNameToIndex[ebn] = 0;
00306 
00307 #ifdef ALBANY_SEACAS
00308   stk::io::put_io_part_attribute(*partVec[0]);
00309 #endif
00310 
00311 
00312   std::vector<std::string> nsNames;
00313   std::string nsn="Bottom";
00314   nsNames.push_back(nsn);
00315   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00316 #ifdef ALBANY_SEACAS
00317     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00318 #endif
00319   nsn="NodeSet0";
00320   nsNames.push_back(nsn);
00321   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00322 #ifdef ALBANY_SEACAS
00323     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00324 #endif
00325   nsn="NodeSet1";
00326   nsNames.push_back(nsn);
00327   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00328 #ifdef ALBANY_SEACAS
00329     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00330 #endif
00331   nsn="NodeSet2";
00332   nsNames.push_back(nsn);
00333   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00334 #ifdef ALBANY_SEACAS
00335     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00336 #endif
00337   nsn="NodeSet3";
00338   nsNames.push_back(nsn);
00339   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00340 #ifdef ALBANY_SEACAS
00341     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00342 #endif
00343   nsn="NodeSet4";
00344   nsNames.push_back(nsn);
00345   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00346 #ifdef ALBANY_SEACAS
00347     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00348 #endif
00349   nsn="NodeSet5";
00350   nsNames.push_back(nsn);
00351   nsPartVec[nsn] = & metaData->declare_part(nsn, metaData->node_rank() );
00352 #ifdef ALBANY_SEACAS
00353     stk::io::put_io_part_attribute(*nsPartVec[nsn]);
00354 #endif
00355 
00356 
00357   std::vector<std::string> ssNames;
00358   std::string ssn="Basal";
00359   ssNames.push_back(ssn);
00360     ssPartVec[ssn] = & metaData->declare_part(ssn, metaData->side_rank() );
00361 #ifdef ALBANY_SEACAS
00362     stk::io::put_io_part_attribute(*ssPartVec[ssn]);
00363 #endif
00364 
00365   stk::mesh::fem::set_cell_topology<shards::Hexahedron<8> >(*partVec[0]);
00366   stk::mesh::fem::set_cell_topology<shards::Quadrilateral<4> >(*ssPartVec[ssn]);
00367 
00368   numDim = 3;
00369   int cub = params->get("Cubature Degree",3);
00370   int worksetSizeMax = params->get("Workset Size",50);
00371   int worksetSize = this->computeWorksetSize(worksetSizeMax, elem_map->NumMyElements());
00372 
00373   const CellTopologyData& ctd = *metaData->get_cell_topology(*partVec[0]).getCellTopologyData();
00374 
00375   this->meshSpecs[0] = Teuchos::rcp(new Albany::MeshSpecsStruct(ctd, numDim, cub,
00376                              nsNames, ssNames, worksetSize, partVec[0]->name(),
00377                              ebNameToIndex, this->interleavedOrdering));
00378 
00379 
00380 }
00381 
00382 Albany::AsciiSTKMeshStruct::~AsciiSTKMeshStruct()
00383 {
00384   delete [] xyz; 
00385   if (have_sh) delete [] sh; 
00386   if (have_bf) delete [] bf; 
00387   delete [] eles; 
00388   delete [] globalElesID; 
00389   delete [] globalNodesID;
00390   delete [] basalFacesID; 
00391 }
00392 
00393 void
00394 Albany::AsciiSTKMeshStruct::setFieldAndBulkData(
00395                                                const Teuchos::RCP<const Epetra_Comm>& comm,
00396                                                const Teuchos::RCP<Teuchos::ParameterList>& params,
00397                                                const unsigned int neq_,
00398                                                const AbstractFieldContainer::FieldContainerRequirements& req,
00399                                                const Teuchos::RCP<Albany::StateInfoStruct>& sis,
00400                                                const unsigned int worksetSize)
00401 {
00402   this->SetupFieldData(comm, neq_, req, sis, worksetSize);
00403 
00404   metaData->commit();
00405 
00406   bulkData->modification_begin(); // Begin modifying the mesh
00407 
00408   stk::mesh::PartVector nodePartVec;
00409   stk::mesh::PartVector singlePartVec(1);
00410   stk::mesh::PartVector emptyPartVec;
00411 #ifdef OUTPUT_TO_SCREEN
00412   *out << "elem_map # elements: " << elem_map->NumMyElements() << std::endl; 
00413   *out << "node_map # elements: " << node_map->NumMyElements() << std::endl;
00414 #endif 
00415   unsigned int ebNo = 0; //element block #??? 
00416   int sideID = 0;
00417 
00418   AbstractSTKFieldContainer::VectorFieldType* coordinates_field = fieldContainer->getCoordinatesField();
00419   AbstractSTKFieldContainer::ScalarFieldType* surfaceHeight_field = fieldContainer->getSurfaceHeightField();
00420   AbstractSTKFieldContainer::ScalarFieldType* flowFactor_field = fieldContainer->getFlowFactorField();
00421   AbstractSTKFieldContainer::ScalarFieldType* temperature_field = fieldContainer->getTemperatureField();
00422   AbstractSTKFieldContainer::ScalarFieldType* basal_friction_field = fieldContainer->getBasalFrictionField();
00423 
00424   if(!surfaceHeight_field) 
00425      have_sh = false;
00426   if(!flowFactor_field) 
00427      have_flwa = false;
00428   if(!temperature_field) 
00429      have_temp = false;
00430   if(!basal_friction_field) 
00431      have_beta = false;
00432 
00433   for (int i=0; i<elem_map->NumMyElements(); i++) {
00434      const unsigned int elem_GID = elem_map->GID(i);
00435      //std::cout << "elem_GID: " << elem_GID << std::endl; 
00436      stk::mesh::EntityId elem_id = (stk::mesh::EntityId) elem_GID;
00437      singlePartVec[0] = partVec[ebNo];
00438      stk::mesh::Entity& elem  = bulkData->declare_entity(metaData->element_rank(), 1+elem_id, singlePartVec);
00439      //I am assuming the ASCII mesh is 1-based not 0-based, so no need to add 1 for STK mesh 
00440      stk::mesh::Entity& llnode = bulkData->declare_entity(metaData->node_rank(), eles[i][0], nodePartVec);
00441      stk::mesh::Entity& lrnode = bulkData->declare_entity(metaData->node_rank(), eles[i][1], nodePartVec);
00442      stk::mesh::Entity& urnode = bulkData->declare_entity(metaData->node_rank(), eles[i][2], nodePartVec);
00443      stk::mesh::Entity& ulnode = bulkData->declare_entity(metaData->node_rank(), eles[i][3], nodePartVec);
00444      stk::mesh::Entity& llnodeb = bulkData->declare_entity(metaData->node_rank(), eles[i][4], nodePartVec);
00445      stk::mesh::Entity& lrnodeb = bulkData->declare_entity(metaData->node_rank(), eles[i][5], nodePartVec);
00446      stk::mesh::Entity& urnodeb = bulkData->declare_entity(metaData->node_rank(), eles[i][6], nodePartVec);
00447      stk::mesh::Entity& ulnodeb = bulkData->declare_entity(metaData->node_rank(), eles[i][7], nodePartVec);
00448      bulkData->declare_relation(elem, llnode, 0);
00449      bulkData->declare_relation(elem, lrnode, 1);
00450      bulkData->declare_relation(elem, urnode, 2);
00451      bulkData->declare_relation(elem, ulnode, 3);
00452      bulkData->declare_relation(elem, llnodeb, 4);
00453      bulkData->declare_relation(elem, lrnodeb, 5);
00454      bulkData->declare_relation(elem, urnodeb, 6);
00455      bulkData->declare_relation(elem, ulnodeb, 7);
00456     
00457 
00458      double* coord;
00459      int node_GID;
00460      unsigned int node_LID;
00461 
00462      node_GID = eles[i][0]-1;
00463      node_LID = node_map->LID(node_GID);
00464      coord = stk::mesh::field_data(*coordinates_field, llnode);
00465      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00466 
00467      node_GID = eles[i][1]-1;
00468      node_LID = node_map->LID(node_GID);
00469      coord = stk::mesh::field_data(*coordinates_field, lrnode);
00470      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00471 
00472      node_GID = eles[i][2]-1;
00473      node_LID = node_map->LID(node_GID);
00474      coord = stk::mesh::field_data(*coordinates_field, urnode);
00475      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00476 
00477      node_GID = eles[i][3]-1;
00478      node_LID = node_map->LID(node_GID);
00479      coord = stk::mesh::field_data(*coordinates_field, ulnode);
00480      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00481 
00482      coord = stk::mesh::field_data(*coordinates_field, llnodeb);
00483      node_GID = eles[i][4]-1;
00484      node_LID = node_map->LID(node_GID);
00485      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00486 
00487      node_GID = eles[i][5]-1;
00488      node_LID = node_map->LID(node_GID);
00489      coord = stk::mesh::field_data(*coordinates_field, lrnodeb);
00490      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00491 
00492      coord = stk::mesh::field_data(*coordinates_field, urnodeb);
00493      node_GID = eles[i][6]-1;
00494      node_LID = node_map->LID(node_GID);
00495      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00496 
00497      coord = stk::mesh::field_data(*coordinates_field, ulnodeb);
00498      node_GID = eles[i][7]-1;
00499      node_LID = node_map->LID(node_GID);
00500      coord[0] = xyz[node_LID][0];   coord[1] = xyz[node_LID][1];   coord[2] = xyz[node_LID][2];
00501 
00502 #ifdef ALBANY_FELIX
00503      if (have_sh) {
00504        double* sHeight;
00505        sHeight = stk::mesh::field_data(*surfaceHeight_field, llnode);
00506        node_GID = eles[i][0]-1;
00507        node_LID = node_map->LID(node_GID);
00508        sHeight[0] = sh[node_LID];
00509 
00510        sHeight = stk::mesh::field_data(*surfaceHeight_field, lrnode);
00511        node_GID = eles[i][1]-1;
00512        node_LID = node_map->LID(node_GID);
00513        sHeight[0] = sh[node_LID];
00514 
00515        sHeight = stk::mesh::field_data(*surfaceHeight_field, urnode);
00516        node_GID = eles[i][2]-1;
00517        node_LID = node_map->LID(node_GID);
00518        sHeight[0] = sh[node_LID];
00519 
00520        sHeight = stk::mesh::field_data(*surfaceHeight_field, ulnode);
00521        node_GID = eles[i][3]-1;
00522        node_LID = node_map->LID(node_GID);
00523        sHeight[0] = sh[node_LID];
00524 
00525        sHeight = stk::mesh::field_data(*surfaceHeight_field, llnodeb);
00526        node_GID = eles[i][4]-1;
00527        node_LID = node_map->LID(node_GID);
00528        sHeight[0] = sh[node_LID];
00529 
00530        sHeight = stk::mesh::field_data(*surfaceHeight_field, lrnodeb);
00531        node_GID = eles[i][5]-1;
00532        node_LID = node_map->LID(node_GID);
00533        sHeight[0] = sh[node_LID];
00534 
00535        sHeight = stk::mesh::field_data(*surfaceHeight_field, urnodeb);
00536        node_GID = eles[i][6]-1;
00537        node_LID = node_map->LID(node_GID);
00538        sHeight[0] = sh[node_LID];
00539 
00540        sHeight = stk::mesh::field_data(*surfaceHeight_field, ulnodeb);
00541        node_GID = eles[i][7]-1;
00542        node_LID = node_map->LID(node_GID);
00543        sHeight[0] = sh[node_LID];
00544      }
00545      if (have_flwa) {
00546        double *flowFactor = stk::mesh::field_data(*flowFactor_field, elem); 
00547        //i is elem_LID (element local ID);
00548        //*out << "i: " << i <<", flwa: " << flwa[i] << std::endl;  
00549        flowFactor[0] = flwa[i]; 
00550      }
00551      if (have_temp) {
00552        double *temperature = stk::mesh::field_data(*temperature_field, elem); 
00553        //i is elem_LID (element local ID);
00554        //*out << "i: " << i <<", temp: " << temperature[i] << std::endl;  
00555        temperature[0] = temper[i]; 
00556      }
00557      if (have_beta) {
00558        double* bFriction; 
00559        bFriction = stk::mesh::field_data(*basal_friction_field, llnode);
00560        node_GID = eles[i][0]-1;
00561        node_LID = node_map->LID(node_GID);
00562        bFriction[0] = beta[node_LID];
00563 
00564        bFriction = stk::mesh::field_data(*basal_friction_field, lrnode);
00565        node_GID = eles[i][1]-1;
00566        node_LID = node_map->LID(node_GID);
00567        bFriction[0] = beta[node_LID];
00568 
00569        bFriction = stk::mesh::field_data(*basal_friction_field, urnode);
00570        node_GID = eles[i][2]-1;
00571        node_LID = node_map->LID(node_GID);
00572        bFriction[0] = beta[node_LID];
00573 
00574        bFriction = stk::mesh::field_data(*basal_friction_field, ulnode);
00575        node_GID = eles[i][3]-1;
00576        node_LID = node_map->LID(node_GID);
00577        bFriction[0] = beta[node_LID];
00578 
00579        bFriction = stk::mesh::field_data(*basal_friction_field, llnodeb);
00580        node_GID = eles[i][4]-1;
00581        node_LID = node_map->LID(node_GID);
00582        bFriction[0] = beta[node_LID];
00583 
00584        bFriction = stk::mesh::field_data(*basal_friction_field, lrnodeb);
00585        node_GID = eles[i][5]-1;
00586        node_LID = node_map->LID(node_GID);
00587        bFriction[0] = beta[node_LID];
00588 
00589        bFriction = stk::mesh::field_data(*basal_friction_field, urnodeb);
00590        node_GID = eles[i][6]-1;
00591        node_LID = node_map->LID(node_GID);
00592        bFriction[0] = beta[node_LID];
00593 
00594        bFriction = stk::mesh::field_data(*basal_friction_field, ulnodeb);
00595        node_GID = eles[i][7]-1;
00596        node_LID = node_map->LID(node_GID);
00597        bFriction[0] = beta[node_LID]; 
00598        }
00599 #endif
00600 
00601      // If first node has z=0 and there is no basal face file provided, identify it as a Basal SS
00602      if (have_bf == false) {
00603        *out <<"No bf file specified...  setting basal boundary to z=0 plane..." << std::endl; 
00604        if ( xyz[eles[i][0]][2] == 0.0) {
00605           //std::cout << "sideID: " << sideID << std::endl; 
00606           singlePartVec[0] = ssPartVec["Basal"];
00607           stk::mesh::EntityId side_id = (stk::mesh::EntityId)(sideID);
00608           sideID++;
00609 
00610          stk::mesh::Entity& side  = bulkData->declare_entity(metaData->side_rank(), 1 + side_id, singlePartVec);
00611          bulkData->declare_relation(elem, side,  4 /*local side id*/);
00612 
00613          bulkData->declare_relation(side, llnode, 0);
00614          bulkData->declare_relation(side, ulnode, 3);
00615          bulkData->declare_relation(side, urnode, 2);
00616          bulkData->declare_relation(side, lrnode, 1);
00617        }
00618     }
00619 
00620     if (xyz[eles[i][0]-1][0] == 0.0) {
00621        singlePartVec[0] = nsPartVec["NodeSet0"];
00622        bulkData->change_entity_parts(llnode, singlePartVec); // node 0
00623        bulkData->change_entity_parts(ulnode, singlePartVec); // 3
00624        bulkData->change_entity_parts(llnodeb, singlePartVec); // 4
00625        bulkData->change_entity_parts(ulnodeb, singlePartVec); // 7
00626     }
00627     if (xyz[eles[i][1]-1][0] == 1.0) {
00628        singlePartVec[0] = nsPartVec["NodeSet1"];
00629        bulkData->change_entity_parts(lrnode, singlePartVec); // 1
00630        bulkData->change_entity_parts(urnode, singlePartVec); // 2
00631        bulkData->change_entity_parts(lrnodeb, singlePartVec); // 5
00632        bulkData->change_entity_parts(urnodeb, singlePartVec); // 6
00633     }
00634     if (xyz[eles[i][0]-1][1] == 0.0) {
00635        singlePartVec[0] = nsPartVec["NodeSet2"];
00636        bulkData->change_entity_parts(llnode, singlePartVec); // 0
00637        bulkData->change_entity_parts(lrnode, singlePartVec); // 1
00638        bulkData->change_entity_parts(llnodeb, singlePartVec); // 4
00639        bulkData->change_entity_parts(lrnodeb, singlePartVec); // 5
00640     }
00641     if (xyz[eles[i][2]-1][1] == 1.0) {
00642        singlePartVec[0] = nsPartVec["NodeSet3"];
00643        bulkData->change_entity_parts(ulnode, singlePartVec); // 3
00644        bulkData->change_entity_parts(urnode, singlePartVec); // 2
00645        bulkData->change_entity_parts(ulnodeb, singlePartVec); // 7
00646        bulkData->change_entity_parts(urnodeb, singlePartVec); // 6
00647     }
00648     if (xyz[eles[i][0]-1][2] == 0.0) {
00649        singlePartVec[0] = nsPartVec["NodeSet4"];
00650        bulkData->change_entity_parts(llnode, singlePartVec); // 0
00651        bulkData->change_entity_parts(lrnode, singlePartVec); // 1
00652        bulkData->change_entity_parts(ulnode, singlePartVec); // 3
00653        bulkData->change_entity_parts(urnode, singlePartVec); // 2
00654     }
00655     if (xyz[eles[i][4]-1][2] == 1.0) {
00656        singlePartVec[0] = nsPartVec["NodeSet5"];
00657        bulkData->change_entity_parts(llnodeb, singlePartVec); // 4
00658        bulkData->change_entity_parts(lrnodeb, singlePartVec); // 5
00659        bulkData->change_entity_parts(ulnodeb, singlePartVec); // 7
00660        bulkData->change_entity_parts(urnodeb, singlePartVec); // 6
00661     }
00662 
00663      // If first node has z=0, identify it as a Bottom NS
00664      if ( xyz[eles[i][0]-1][2] == 0.0) {
00665        singlePartVec[0] = nsPartVec["Bottom"];
00666        bulkData->change_entity_parts(llnode, singlePartVec); // 0
00667        bulkData->change_entity_parts(lrnode, singlePartVec); // 1
00668        bulkData->change_entity_parts(ulnode, singlePartVec); // 3
00669        bulkData->change_entity_parts(urnode, singlePartVec); // 2
00670      }
00671   }
00672   if (have_bf == true) {
00673     *out << "Setting basal surface connectivity from bf file provided..." << std::endl;  
00674     for (int i=0; i<basal_face_map->NumMyElements(); i++) {
00675        singlePartVec[0] = ssPartVec["Basal"];
00676        sideID = basal_face_map->GID(i); 
00677        stk::mesh::EntityId side_id = (stk::mesh::EntityId)(sideID);
00678        stk::mesh::Entity& side  = bulkData->declare_entity(metaData->side_rank(), 1 + side_id, singlePartVec);
00679 
00680        const unsigned int elem_GID = bf[i][0];
00681        //std::cout << "elem_GID: " << elem_GID << std::endl; 
00682        stk::mesh::EntityId elem_id = (stk::mesh::EntityId) elem_GID;
00683        stk::mesh::Entity& elem  = bulkData->declare_entity(metaData->element_rank(), elem_id, emptyPartVec);
00684        bulkData->declare_relation(elem, side,  4 /*local side id*/);
00685        stk::mesh::Entity& llnode = bulkData->declare_entity(metaData->node_rank(), bf[i][1], nodePartVec);
00686        stk::mesh::Entity& lrnode = bulkData->declare_entity(metaData->node_rank(), bf[i][2], nodePartVec);
00687        stk::mesh::Entity& urnode = bulkData->declare_entity(metaData->node_rank(), bf[i][3], nodePartVec);
00688        stk::mesh::Entity& ulnode = bulkData->declare_entity(metaData->node_rank(), bf[i][4], nodePartVec);
00689        
00690        bulkData->declare_relation(side, llnode, 0);
00691        bulkData->declare_relation(side, ulnode, 3);
00692        bulkData->declare_relation(side, urnode, 2);
00693        bulkData->declare_relation(side, lrnode, 1);
00694     }
00695   }
00696 
00697   bulkData->modification_end();
00698 }
00699 
00700 Teuchos::RCP<const Teuchos::ParameterList>
00701 Albany::AsciiSTKMeshStruct::getValidDiscretizationParameters() const
00702 {
00703   Teuchos::RCP<Teuchos::ParameterList> validPL =
00704     this->getValidGenericSTKParameters("Valid ASCII_DiscParams");
00705 
00706   return validPL;
00707 }