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

Go to the documentation of this file.
00001 // ===================================================
00003 // ===================================================
00004 
00005 #include "Interface.hpp"
00006 #include "Albany_MpasSTKMeshStruct.hpp"
00007 #include "Teuchos_ParameterList.hpp"
00008 #include "Teuchos_RCP.hpp"
00009 #include "Albany_Utils.hpp"
00010 #include "Albany_SolverFactory.hpp"
00011 #include "Teuchos_XMLParameterListHelpers.hpp"
00012 #include <stk_mesh/base/FieldData.hpp>
00013 #include "Piro_PerformSolve.hpp"
00014 #include "Thyra_EpetraThyraWrappers.hpp"
00015 #include <stk_io/IossBridge.hpp>
00016 #include <stk_io/MeshReadWriteUtils.hpp>
00017 #include <stk_mesh/base/GetEntities.hpp>
00018 #include <stk_mesh/base/FieldData.hpp>
00019 #include <Ionit_Initializer.h>
00020 #include "Albany_OrdinarySTKFieldContainer.hpp"
00021 
00022 
00023 
00024 
00025 
00026 // ===================================================
00028 // ===================================================
00029 //using namespace LifeV;
00030 
00031 //typedef std::list<exchange> exchangeList_Type;
00032 
00033 // ice_problem pointer
00034 //ICEProblem *iceProblemPtr = 0;
00035 Teuchos::RCP<Albany::MpasSTKMeshStruct> meshStruct2D;
00036 Teuchos::RCP<Albany::MpasSTKMeshStruct> meshStruct;
00037 Teuchos::RCP<const Epetra_Comm> mpiComm;
00038 Teuchos::RCP<Teuchos::ParameterList> appParams;
00039 Teuchos::RCP<Teuchos::ParameterList> discParams;
00040 Teuchos::RCP<Albany::SolverFactory> slvrfctry;
00041 Teuchos::RCP<Thyra::ModelEvaluator<double> > solver;
00042 int Ordering =0; //ordering ==0 means that the mesh is extruded layerwise, whereas ordering==1 means that the mesh is extruded columnwise.
00043 MPI_Comm comm, reducedComm;
00044 bool isDomainEmpty = true;
00045 bool initialize_velocity = true;
00046 bool first_time_step = true;
00047 int nCells_F, nEdges_F, nVertices_F;
00048 int nCellsSolve_F, nEdgesSolve_F, nVerticesSolve_F;
00049 int nVertices, nEdges, nTriangles, nGlobalVertices, nGlobalEdges, nGlobalTriangles;
00050 int maxNEdgesOnCell_F;
00051 int const *cellsOnEdge_F, *cellsOnVertex_F, *verticesOnCell_F, *verticesOnEdge_F, *edgesOnCell_F, *indexToCellID_F, *nEdgesOnCells_F;
00052 std::vector<double> layersRatio, levelsNormalizedThickness;
00053 int nLayers;
00054 double const *xCell_F, *yCell_F, *zCell_F, *areaTriangle_F;
00055 std::vector<double> xCellProjected, yCellProjected, zCellProjected;
00056 const double unit_length = 1000;
00057 const double T0 = 273.15;
00058 const double minThick =1e-2; //10m
00059 const double minBeta =1e-5;
00060 void *phgGrid = 0;
00061 std::vector<int> edgesToReceive, fCellsToReceive, indexToTriangleID, verticesOnTria, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge;
00062 std::vector<int> indexToVertexID, vertexToFCell, indexToEdgeID, edgeToFEdge, mask, fVertexToTriangleID, fCellToVertex;
00063 std::vector<double> temperatureOnTetra, velocityOnVertices, velocityOnCells, elevationData, thicknessData, betaData, smb_F, thicknessOnCells;
00064 std::vector<bool> isVertexBoundary, isBoundaryEdge;
00065 ;
00066 int numBoundaryEdges;
00067 double radius;
00068 
00069 exchangeList_Type const *sendCellsList_F=0, *recvCellsList_F=0;
00070 exchangeList_Type const *sendEdgesList_F=0, *recvEdgesList_F=0;
00071 exchangeList_Type const *sendVerticesList_F=0, *recvVerticesList_F=0;
00072 exchangeList_Type sendCellsListReversed, recvCellsListReversed, sendEdgesListReversed, recvEdgesListReversed;
00073 
00074 typedef struct TET_ {
00075     int verts[4];
00076     int neighbours[4];
00077     char bound_type[4];
00078 } TET;
00079 
00080 
00081 
00082 exchange::exchange(int _procID, int const *  vec_first, int const *  vec_last, int fieldDim):
00083         procID(_procID),
00084         vec(vec_first, vec_last),
00085         buffer(fieldDim*(vec_last-vec_first)),
00086         doubleBuffer(fieldDim*(vec_last-vec_first)){}
00087 
00088 
00089 /***********************************************************/
00090 // epetra <-> thyra conversion utilities
00091 Teuchos::RCP<const Epetra_Vector>
00092 epetraVectorFromThyra(
00093   const Teuchos::RCP<const Epetra_Comm> &comm,
00094   const Teuchos::RCP<const Thyra::VectorBase<double> > &thyra)
00095 {
00096   Teuchos::RCP<const Epetra_Vector> result;
00097   if (Teuchos::nonnull(thyra)) {
00098     const Teuchos::RCP<const Epetra_Map> epetra_map = Thyra::get_Epetra_Map(*thyra->space(), comm);
00099     result = Thyra::get_Epetra_Vector(*epetra_map, thyra);
00100   }
00101   return result;
00102 }
00103 
00104 Teuchos::RCP<const Epetra_MultiVector>
00105 epetraMultiVectorFromThyra(
00106   const Teuchos::RCP<const Epetra_Comm> &comm,
00107   const Teuchos::RCP<const Thyra::MultiVectorBase<double> > &thyra)
00108 {
00109   Teuchos::RCP<const Epetra_MultiVector> result;
00110   if (Teuchos::nonnull(thyra)) {
00111     const Teuchos::RCP<const Epetra_Map> epetra_map = Thyra::get_Epetra_Map(*thyra->range(), comm);
00112     result = Thyra::get_Epetra_MultiVector(*epetra_map, thyra);
00113   }
00114   return result;
00115 }
00116 
00117 void epetraFromThyra(
00118   const Teuchos::RCP<const Epetra_Comm> &comm,
00119   const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > &thyraResponses,
00120   const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > &thyraSensitivities,
00121   Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > &responses,
00122   Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > &sensitivities)
00123 {
00124   responses.clear();
00125   responses.reserve(thyraResponses.size());
00126   typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > ThyraResponseArray;
00127   for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
00128       it_end = thyraResponses.end(),
00129       it = it_begin;
00130       it != it_end;
00131       ++it) {
00132     responses.push_back(epetraVectorFromThyra(comm, *it));
00133   }
00134 
00135   sensitivities.clear();
00136   sensitivities.reserve(thyraSensitivities.size());
00137   typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > ThyraSensitivityArray;
00138   for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
00139       it_end = thyraSensitivities.end(),
00140       it = it_begin;
00141       it != it_end;
00142       ++it) {
00143     ThyraSensitivityArray::const_reference sens_thyra = *it;
00144     Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > sens;
00145     sens.reserve(sens_thyra.size());
00146     for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
00147         jt_end = sens_thyra.end();
00148         jt != jt_end;
00149         ++jt) {
00150         sens.push_back(epetraMultiVectorFromThyra(comm, *jt));
00151     }
00152     sensitivities.push_back(sens);
00153   }
00154 }
00155 
00156 /***********************************************************/
00157 
00158 extern "C"
00159 {
00160 
00161 #ifdef HAVE_PHG
00162     extern void phg_init_(int *fComm); 
00163     extern void *phgImportParallelGrid(void *old_grid,
00164               int nvert,
00165               int nelem, 
00166               int nvert_global,
00167               int nelem_global,
00168               int *L2Gmap_vert,
00169               int *L2Gmap_elem,
00170               double *coord,
00171               TET *tet,
00172               MPI_Comm comm
00173               );
00174     extern void phgSolveIceStokes(void *g_ptr, 
00175           int nlayer,
00176           const double *T,
00177           const double *beta,
00178           double *U);
00179 #endif
00180 
00181   // ===================================================
00183   // ===================================================
00184     int velocity_solver_init_mpi(int *fComm)
00185     {
00186 
00187     // get MPI_Comm from Fortran
00188     comm = MPI_Comm_f2c(*fComm);
00189 
00190     return 0;
00191   }
00192 
00193 
00194   void velocity_solver_set_grid_data(int const * _nCells_F, int const * _nEdges_F, int const * _nVertices_F, int const * _nLayers,
00195                                  int const * _nCellsSolve_F, int const * _nEdgesSolve_F, int const * _nVerticesSolve_F, int const* _maxNEdgesOnCell_F,
00196                                  double const * radius_F,
00197                                  int const * _cellsOnEdge_F, int const * _cellsOnVertex_F, int const * _verticesOnCell_F, int const * _verticesOnEdge_F, int const * _edgesOnCell_F,
00198                                  int const* _nEdgesOnCells_F, int const * _indexToCellID_F,
00199                                  double const *  _xCell_F, double const *  _yCell_F, double const *  _zCell_F, double const *  _areaTriangle_F,
00200                                  int const * sendCellsArray_F, int const * recvCellsArray_F,
00201                                  int const * sendEdgesArray_F, int const * recvEdgesArray_F,
00202                                  int const * sendVerticesArray_F, int const * recvVerticesArray_F)
00203   {
00204       nCells_F = *_nCells_F;
00205       nEdges_F = *_nEdges_F;
00206       nVertices_F = *_nVertices_F;
00207       nLayers = *_nLayers;
00208       nCellsSolve_F = *_nCellsSolve_F;
00209       nEdgesSolve_F = *_nEdgesSolve_F;
00210       nVerticesSolve_F = *_nVerticesSolve_F;
00211       maxNEdgesOnCell_F = *_maxNEdgesOnCell_F;
00212       radius = *radius_F;
00213       cellsOnEdge_F = _cellsOnEdge_F;
00214       cellsOnVertex_F = _cellsOnVertex_F;
00215       verticesOnCell_F = _verticesOnCell_F;
00216       verticesOnEdge_F = _verticesOnEdge_F;
00217       edgesOnCell_F = _edgesOnCell_F;
00218       nEdgesOnCells_F =_nEdgesOnCells_F;
00219       indexToCellID_F = _indexToCellID_F;
00220       xCell_F = _xCell_F;
00221       yCell_F = _yCell_F;
00222       zCell_F = _zCell_F;
00223       areaTriangle_F = _areaTriangle_F;
00224       mask.resize(nVertices_F);
00225 
00226       thicknessOnCells.resize(nCellsSolve_F);
00227 
00228       sendCellsList_F = new exchangeList_Type(unpackMpiArray(sendCellsArray_F));
00229       recvCellsList_F = new exchangeList_Type(unpackMpiArray(recvCellsArray_F));
00230       sendEdgesList_F = new exchangeList_Type(unpackMpiArray(sendEdgesArray_F));
00231       recvEdgesList_F = new exchangeList_Type(unpackMpiArray(recvEdgesArray_F));
00232       sendVerticesList_F = new exchangeList_Type(unpackMpiArray(sendVerticesArray_F));
00233       recvVerticesList_F = new exchangeList_Type(unpackMpiArray(recvVerticesArray_F));
00234 
00235       if(radius > 10)
00236       {
00237         xCellProjected.resize(nCells_F);
00238         yCellProjected.resize(nCells_F);
00239         zCellProjected.assign(nCells_F, 0.);
00240         for(int i=0; i<nCells_F; i++)
00241         {
00242           double r=std::sqrt(xCell_F[i]*xCell_F[i]+yCell_F[i]*yCell_F[i]+zCell_F[i]*zCell_F[i]);
00243           xCellProjected[i] = radius*std::asin(xCell_F[i]/r);
00244           yCellProjected[i] = radius*std::asin(yCell_F[i]/r);
00245         }
00246         xCell_F=&xCellProjected[0];
00247         yCell_F=&yCellProjected[0];
00248         zCell_F=&zCellProjected[0];
00249       }
00250   }
00251 
00252   void velocity_solver_init_l1l2(double const * levelsRatio_F)
00253   {
00254       velocityOnVertices.resize(2*nVertices*(nLayers+1),0.);
00255       velocityOnCells.resize(2*nCells_F*(nLayers+1),0.);
00256 
00257     if(isDomainEmpty)
00258         return;
00259 
00260     layersRatio.resize(nLayers);
00261     // !!Indexing of layers is reversed
00262     for(int i=0; i<nLayers; i++)
00263         layersRatio[i] = levelsRatio_F[nLayers-1-i];
00264     //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin());
00265     mapCellsToVertices(velocityOnCells, velocityOnVertices, 2, nLayers, Ordering);
00266     //std::cout << __LINE__ << " max: " << *std::max_element(velocityOnVertices.begin(), velocityOnVertices.end()) <<", min: " << *std::min_element(velocityOnVertices.begin(), velocityOnVertices.end()) <<std::endl;
00267    //   iceProblemPtr->initializeSolverL1L2(layersRatio, velocityOnVertices, initialize_velocity);
00268       initialize_velocity = false;
00269   }
00270 
00271 
00272 
00273 #define GET_HERE printf("get %25s, line: %d\n", __FILE__, __LINE__);
00274 
00275 
00276 
00277 
00278 
00279 void velocity_solver_export_2d_data(double const * lowerSurface_F, double const * thickness_F,
00280                             double const * beta_F)
00281 {
00282       if(isDomainEmpty)
00283       return;
00284 
00285       import2DFields(lowerSurface_F, thickness_F, beta_F, minThick);
00286 
00287       Teuchos::RCP<stk::io::MeshData> mesh_data =Teuchos::rcp(new stk::io::MeshData);
00288       stk::io::create_output_mesh("mesh2D.exo", reducedComm, *meshStruct2D->bulkData, *mesh_data);
00289       stk::io::define_output_fields(*mesh_data, *meshStruct->metaData);
00290 
00291       //  iceProblemPtr->export_2D_fields(elevationData, thicknessData, betaData, indexToVertexID);
00292 }
00293 
00294 
00295 
00296 void velocity_solver_solve_l1l2(double const * lowerSurface_F, double const * thickness_F,
00297                             double const * beta_F, double const * temperature_F,
00298                             double * u_normal_F,
00299                             double * /*heatIntegral_F*/, double * /*viscosity_F*/)
00300   {
00301   }
00302 
00303 
00304   void velocity_solver_export_l1l2_velocity()
00305     {
00306       if(isDomainEmpty)
00307           return;
00308 
00309    /*   if(iceProblemPtr->mesh3DPtr == 0)
00310       {
00311           levelsNormalizedThickness.resize(nLayers+1);
00312 
00313           levelsNormalizedThickness[0] = 0;
00314           for(int i=0; i< nLayers; i++)
00315               levelsNormalizedThickness[i+1] = levelsNormalizedThickness[i]+layersRatio[i];
00316 
00317 
00318           iceProblemPtr->mesh3DPtr.reset(new RegionMesh<LinearTetra>() );
00319 
00320           std::vector<double> regulThk(thicknessData);
00321           for(int index=0; index<nVertices; index++)
00322             regulThk[index] = std::max(1e-4,thicknessData[index]);
00323 
00324           std::vector<int> mpasIndexToVertexID(nVertices);
00325           for(int i=0; i<nVertices; i++)
00326             mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]];
00327 
00328           meshVerticallyStructuredFrom2DMesh( *(iceProblemPtr->mesh3DPtr), *(iceProblemPtr->mesh2DPtr), mpasIndexToVertexID, levelsNormalizedThickness,
00329                    elevationData, regulThk, LayerWise, 1, 2, reducedComm );
00330       }
00331         iceProblemPtr->export_solution_L1L2();
00332 
00333         */
00334     }
00335 
00336 
00337   void velocity_solver_init_fo(double const * levelsRatio_F)
00338     {
00339       velocityOnVertices.resize(2*nVertices*(nLayers+1),0.);
00340       velocityOnCells.resize(2*nCells_F*(nLayers+1),0.);
00341 
00342         if(isDomainEmpty)
00343             return;
00344 
00345         layersRatio.resize(nLayers);
00346         // !!Indexing of layers is reversed
00347         for(int i=0; i<nLayers; i++)
00348             layersRatio[i] = levelsRatio_F[nLayers-1-i];
00349         //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin());
00350 
00351         mapCellsToVertices(velocityOnCells, velocityOnVertices, 2, nLayers, Ordering);
00352     //    iceProblemPtr->initializeSolverFO(layersRatio, velocityOnVertices, thicknessData, elevationData, indexToVertexID, initialize_velocity);
00353         initialize_velocity = false;
00354     }
00355 
00356 
00357   void velocity_solver_solve_fo(double const * lowerSurface_F, double const * thickness_F,
00358                                 double const * beta_F, double const * temperature_F,
00359                                 double * u_normal_F,
00360                                 double * /*heatIntegral_F*/, double * /*viscosity_F*/)
00361     {
00362     std::fill (u_normal_F, u_normal_F + nEdges_F * nLayers, 0.);
00363 
00364     if(!isDomainEmpty)
00365     {
00366 
00367           import2DFields(lowerSurface_F, thickness_F, beta_F, minThick);
00368           importP0Temperature(temperature_F);
00369 
00370           std::vector<double> regulThk(thicknessData);
00371       for(int index=0; index<nVertices; index++)
00372         regulThk[index] = std::max(1e-4,thicknessData[index]);
00373 
00374 
00375       UInt numVertices3D =  (nLayers+1)*indexToVertexID.size();
00376       UInt numPrisms =  nLayers*indexToTriangleID.size();
00377       int vertexColumnShift = (Ordering == 1) ? 1 : nGlobalVertices;
00378         int lVertexColumnShift = (Ordering == 1) ? 1 : indexToVertexID.size();
00379       int vertexLayerShift = (Ordering == 0) ? 1 : nLayers+1;
00380 
00381       int elemColumnShift = (Ordering == 1) ? 3 : 3*nGlobalTriangles;
00382       int lElemColumnShift = (Ordering == 1) ? 3 : 3*indexToTriangleID.size();
00383       int elemLayerShift = (Ordering == 0) ? 3 : 3*nLayers;
00384 
00385       const bool interleavedOrdering = meshStruct->getInterleavedOrdering();
00386       Albany::AbstractSTKFieldContainer::VectorFieldType* solutionField;
00387 
00388       if(interleavedOrdering)
00389         solutionField = Teuchos::rcp_dynamic_cast<Albany::OrdinarySTKFieldContainer<true> >(meshStruct->getFieldContainer())->getSolutionField();
00390       else
00391         solutionField = Teuchos::rcp_dynamic_cast<Albany::OrdinarySTKFieldContainer<false> >(meshStruct->getFieldContainer())->getSolutionField();
00392 
00393       for ( UInt j = 0 ; j < numVertices3D ; ++j )
00394         {
00395          int ib = (Ordering == 0)*(j%lVertexColumnShift) + (Ordering == 1)*(j/vertexLayerShift);
00396          int il = (Ordering == 0)*(j/lVertexColumnShift) + (Ordering == 1)*(j%vertexLayerShift);
00397          int gId = il*vertexColumnShift+vertexLayerShift * indexToVertexID[ib];
00398          stk::mesh::Entity& node = *meshStruct->bulkData->get_entity(meshStruct->metaData->node_rank(), gId+1);
00399          double* coord = stk::mesh::field_data(*meshStruct->getCoordinatesField(), node);
00400          coord[2] = elevationData[ib] - levelsNormalizedThickness[nLayers-il]*regulThk[ib];
00401          double* sHeight = stk::mesh::field_data(*meshStruct->getFieldContainer()->getSurfaceHeightField(), node);
00402          sHeight[0] = elevationData[ib];
00403          double* thickness = stk::mesh::field_data(*meshStruct->getFieldContainer()->getThicknessField(),node);
00404          thickness[0] = thicknessData[ib];
00405          double* sol = stk::mesh::field_data(*solutionField, node);
00406          sol[0] = velocityOnVertices[j];
00407          sol[1] = velocityOnVertices[j + numVertices3D];
00408          if(il ==0)
00409          {
00410            double* beta = stk::mesh::field_data(*meshStruct->getFieldContainer()->getBasalFrictionField(),node);
00411            beta[0] = std::max(betaData[ib], minBeta);
00412          }
00413         }
00414 
00415 
00416       for ( UInt j = 0 ; j < numPrisms ; ++j )
00417       {
00418          int ib = (Ordering == 0)*(j%(lElemColumnShift/3)) + (Ordering == 1)*(j/(elemLayerShift/3));
00419          int il = (Ordering == 0)*(j/(lElemColumnShift/3)) + (Ordering == 1)*(j%(elemLayerShift/3));
00420          int gId = il*elemColumnShift+elemLayerShift * indexToTriangleID[ib];
00421          int lId = il*lElemColumnShift+elemLayerShift * ib;
00422          for(int iTetra=0; iTetra<3; iTetra++)
00423          {
00424            stk::mesh::Entity& elem = *meshStruct->bulkData->get_entity(meshStruct->metaData->element_rank(), ++gId);
00425            double* temperature = stk::mesh::field_data(*meshStruct->getFieldContainer()->getTemperatureField(), elem);
00426            temperature[0] = temperatureOnTetra[lId++] ;
00427          }
00428       }
00429 
00430       meshStruct->setHasRestartSolution(!first_time_step);
00431 
00432       Teuchos::RCP<Albany::AbstractSTKMeshStruct> stkMeshStruct = meshStruct;
00433       discParams->set("STKMeshStruct",stkMeshStruct);
00434       Teuchos::RCP<Teuchos::ParameterList> paramList = Teuchos::rcp(&slvrfctry->getParameters(),false);
00435       if(!first_time_step)
00436       {
00437         meshStruct->setRestartDataTime(paramList->sublist("Problem").get("Homotopy Restart Step", 1.));
00438         double homotopy = paramList->sublist("Problem").sublist("FELIX Viscosity").get("Glen's Law Homotopy Parameter", 1.0);
00439         if(meshStruct->restartDataTime()== homotopy)
00440           paramList->sublist("Problem").set("Solution Method", "Steady");
00441       }
00442       Teuchos::RCP<Albany::Application> app = Teuchos::rcp(new Albany::Application(mpiComm, paramList));
00443       solver = slvrfctry->createThyraSolverAndGetAlbanyApp(app, mpiComm, mpiComm);
00444 
00445 
00446        Teuchos::ParameterList solveParams;
00447        solveParams.set("Compute Sensitivities", false);
00448 
00449        Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > thyraResponses;
00450        Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > thyraSensitivities;
00451        Piro::PerformSolveBase(*solver, solveParams, thyraResponses, thyraSensitivities);
00452 
00453        //  Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > responses;
00454        //  Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > sensitivities;
00455        //  epetraFromThyra(mpiComm, thyraResponses, thyraSensitivities, responses, sensitivities);
00456 
00457        // get solution vector out
00458        //const Teuchos::RCP<const Epetra_Vector> solution = responses.back();
00459 
00460        const Epetra_Map& overlapMap(*app->getDiscretization()->getOverlapMap());
00461        Epetra_Import import(overlapMap, *app->getDiscretization()->getMap());
00462        Epetra_Vector solution(overlapMap);
00463        solution.Import(*app->getDiscretization()->getSolutionField(), import, Insert);
00464 
00465        //UInt componentGlobalLength = (nLayers+1)*nGlobalVertices; //mesh3DPtr->numGlobalVertices();
00466 
00467        for ( UInt j = 0 ; j < numVertices3D ; ++j )
00468        {
00469          int ib = (Ordering == 0)*(j%lVertexColumnShift) + (Ordering == 1)*(j/vertexLayerShift);
00470          int il = (Ordering == 0)*(j/lVertexColumnShift) + (Ordering == 1)*(j%vertexLayerShift);
00471          int gId = il*vertexColumnShift+vertexLayerShift * indexToVertexID[ib];
00472         
00473          int lId0, lId1;
00474 
00475          if(interleavedOrdering)
00476          {
00477            lId0= overlapMap.LID(2*gId);
00478            lId1 = lId0+1;
00479          }
00480          else
00481          {
00482            lId0 = overlapMap.LID(gId);
00483            lId1 = lId0+numVertices3D;
00484          }
00485          velocityOnVertices[j] = solution[lId0];
00486          velocityOnVertices[j + numVertices3D] = solution[lId1];
00487        }
00488 
00489        std::vector<int> mpasIndexToVertexID (nVertices);
00490        for (int i = 0; i < nVertices; i++)
00491        {
00492          mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]];
00493        }
00494        get_tetraP1_velocity_on_FEdges (u_normal_F, velocityOnVertices, edgeToFEdge, mpasIndexToVertexID);
00495     }
00496 
00497 
00498     mapVerticesToCells (velocityOnVertices, &velocityOnCells[0], 2, nLayers, Ordering);
00499 
00500 
00501     std::vector<double> velOnEdges (nEdges * nLayers);
00502     for (int i = 0; i < nEdges; i++)
00503     {
00504       for (int il = 0; il < nLayers; il++)
00505       {
00506         velOnEdges[i * nLayers + il] = u_normal_F[edgeToFEdge[i] * nLayers + il];
00507       }
00508     }
00509 
00510     allToAll (u_normal_F,  &sendEdgesListReversed, &recvEdgesListReversed, nLayers);
00511 
00512     allToAll (u_normal_F,  sendEdgesList_F, recvEdgesList_F, nLayers);
00513 
00514     first_time_step = false;
00515    }
00516 
00517 
00518 
00519   void velocity_solver_export_fo_velocity()
00520   {
00521     if(isDomainEmpty)
00522       return;
00523 
00524     Ioss::Init::Initializer io;
00525       Teuchos::RCP<stk::io::MeshData> mesh_data =Teuchos::rcp(new stk::io::MeshData);
00526       stk::io::create_output_mesh("IceSheet.exo", reducedComm, *meshStruct->bulkData, *mesh_data);
00527       stk::io::define_output_fields(*mesh_data, *meshStruct->metaData);
00528       stk::io::process_output_request(*mesh_data, *meshStruct->bulkData, 0.0);
00529   }
00530 
00531 
00532   void velocity_solver_finalize()
00533   {
00534      // delete iceProblemPtr;
00535       delete sendCellsList_F;
00536         delete recvCellsList_F;
00537         delete sendEdgesList_F;
00538         delete recvEdgesList_F;
00539         delete sendVerticesList_F;
00540         delete recvVerticesList_F;
00541   }
00542 
00543 
00544 
00545   /*duality:
00546    *
00547    *   mpas(F) |  lifev
00548    *  ---------|---------
00549    *   cell    |  vertex
00550    *   vertex  |  triangle
00551    *   edge    |  edge
00552    *
00553    */
00554 
00555   void velocity_solver_compute_2d_grid(int const * verticesMask_F)
00556   {
00557       int numProcs,me ;
00558 
00559 
00560         MPI_Comm_size ( comm, &numProcs );
00561         MPI_Comm_rank ( comm, &me );
00562         std::vector<int> partialOffset(numProcs+1), globalOffsetTriangles(numProcs+1),
00563             globalOffsetVertices(numProcs+1), globalOffsetEdge(numProcs+1);
00564 
00565 
00566         std::vector<int> triangleToFVertex; triangleToFVertex.reserve(nVertices_F);
00567         std::vector<int> fVertexToTriangle(nVertices_F, NotAnId);
00568         bool changed=false;
00569         for(int i(0); i<nVerticesSolve_F; i++)
00570         {
00571             if((verticesMask_F[i] & 0x02)&& !isGhostTriangle(i))
00572             {
00573                 fVertexToTriangle[i] = triangleToFVertex.size();
00574                 triangleToFVertex.push_back(i);
00575             }
00576             changed = changed || (verticesMask_F[i]!=mask[i]);
00577         }
00578 
00579         for(int i(0); i<nVertices_F; i++)
00580           mask[i] = verticesMask_F[i];
00581 
00582         if(changed) std::cout<< "mask changed!!" <<std::endl;
00583 
00584         if((me==0)&&(triangleToFVertex.size()==0))
00585         for(int i(0); i<nVerticesSolve_F; i++)
00586         {
00587             if (!isGhostTriangle(i))
00588             {
00589                 fVertexToTriangle[i] = triangleToFVertex.size();
00590                 triangleToFVertex.push_back(i);
00591                 break;
00592             }
00593         }
00594 
00595         nTriangles = triangleToFVertex.size();
00596 
00597         initialize_iceProblem(nTriangles);
00598 
00599         //Compute the global number of triangles, and the localOffset on the local processor, such that a globalID = localOffset + index
00600         int localOffset(0);
00601         nGlobalTriangles=0;
00602         computeLocalOffset(nTriangles, localOffset, nGlobalTriangles);
00603 
00604 
00605         //Communicate the globalIDs, computed locally, to the other processors.
00606         indexToTriangleID.resize(nTriangles);
00607 
00608         //To make local, not used
00609         fVertexToTriangleID.assign(nVertices_F, NotAnId);
00610      //   std::vector<int> fVertexToTriangleID(nVertices_F, NotAnId);
00611         for(int index(0); index< nTriangles; index++)
00612             fVertexToTriangleID[ triangleToFVertex[index] ] = index+localOffset;
00613 
00614         allToAll(fVertexToTriangleID, sendVerticesList_F, recvVerticesList_F);
00615 
00616         for(int index(0); index< nTriangles; index++)
00617             indexToTriangleID[index] = fVertexToTriangleID[ triangleToFVertex[index] ];
00618 
00619 
00620         //Compute triangle edges
00621         std::vector<int> fEdgeToEdge(nEdges_F), edgesToSend, trianglesProcIds(nVertices_F);
00622         getProcIds(trianglesProcIds, recvVerticesList_F);
00623 
00624         int interfaceSize(0);
00625 
00626         std::vector<int> fEdgeToEdgeID(nEdges_F, NotAnId);
00627         edgesToReceive.clear();
00628     edgeToFEdge.clear();
00629     isBoundaryEdge.clear();
00630     trianglesOnEdge.clear();
00631     
00632 
00633         edgesToReceive.reserve(nEdges_F-nEdgesSolve_F);
00634         edgeToFEdge.reserve(nEdges_F);
00635         trianglesOnEdge.reserve(nEdges_F*2);
00636         edgesToSend.reserve(nEdgesSolve_F);
00637         isBoundaryEdge.reserve(nEdges_F);
00638 
00639 
00640 
00641         //first, we compute boundary edges (boundary edges must be the first edges)
00642         for( int i=0; i<nEdges_F; i++)
00643         {
00644             ID fVertex1(verticesOnEdge_F[2*i]-1), fVertex2(verticesOnEdge_F[2*i+1]-1);
00645             ID triaId_1 = fVertexToTriangleID[fVertex1];
00646             ID triaId_2 = fVertexToTriangleID[fVertex2];
00647             bool isboundary = (triaId_1 == NotAnId) || (triaId_2 == NotAnId);
00648 
00649             ID iTria1 = fVertexToTriangle[fVertex1];
00650             ID iTria2 = fVertexToTriangle[fVertex2];
00651             if(iTria1 == NotAnId) std::swap(iTria1,iTria2);
00652             bool belongsToLocalTriangle = (iTria1 != NotAnId) || (iTria2 != NotAnId);
00653 
00654             if( belongsToLocalTriangle)
00655             {
00656                 if(isboundary )
00657                 {
00658                     fEdgeToEdge[i] = edgeToFEdge.size();
00659                     edgeToFEdge.push_back(i);
00660                     trianglesOnEdge.push_back(iTria1);
00661                     trianglesOnEdge.push_back(iTria2);
00662                     isBoundaryEdge.push_back( true );
00663                 }
00664                 else
00665                     interfaceSize += (iTria2 == NotAnId);
00666             }
00667         }
00668 
00669         numBoundaryEdges = edgeToFEdge.size();
00670 
00671         //procOnInterfaceEdge contains the pairs <id of interface edge, rank of processor that owns the non local element adjacent to the edge>.
00672         std::vector<std::pair<int,int> >  procOnInterfaceEdge;
00673         procOnInterfaceEdge.reserve(interfaceSize);
00674 
00675         //then, we compute the other edges
00676         for(int i=0; i<nEdges_F; i++)
00677         {
00678 
00679             ID fVertex1(verticesOnEdge_F[2*i]-1), fVertex2(verticesOnEdge_F[2*i+1]-1);
00680             ID iTria1 = fVertexToTriangle[fVertex1];
00681             ID iTria2 = fVertexToTriangle[fVertex2];
00682 
00683             ID triaId_1 = fVertexToTriangleID[fVertex1]; //global Triangle
00684             ID triaId_2 = fVertexToTriangleID[fVertex2]; //global Triangle
00685 
00686             if(iTria1 == NotAnId)
00687             {
00688                 std::swap(iTria1,iTria2);
00689                 std::swap(fVertex1,fVertex2);
00690             }
00691 
00692 
00693             bool belongsToAnyTriangle = (triaId_1 != NotAnId) || (triaId_2 != NotAnId);
00694             bool isboundary = (triaId_1 == NotAnId) || (triaId_2 == NotAnId);
00695             bool belongsToLocalTriangle = (iTria1 != NotAnId);
00696             bool isMine = i<nEdgesSolve_F;
00697 
00698             if( belongsToLocalTriangle && !isboundary)
00699             {
00700                 fEdgeToEdge[i] = edgeToFEdge.size();
00701                 edgeToFEdge.push_back(i);
00702                 trianglesOnEdge.push_back(iTria1);
00703                 trianglesOnEdge.push_back(iTria2);
00704                 isBoundaryEdge.push_back( false );
00705                 if(iTria2 == NotAnId)
00706                     procOnInterfaceEdge.push_back(std::make_pair(fEdgeToEdge[i],trianglesProcIds[fVertex2]));
00707             }
00708 
00709             if( belongsToAnyTriangle &&  isMine)
00710             {
00711                 edgesToSend.push_back(i);
00712                 if(! belongsToLocalTriangle)
00713                     edgesToReceive.push_back(i);
00714             }
00715 
00716         }
00717 
00718         //Compute the global number of edges, and the localOffset on the local processor, such that a globalID = localOffset + index
00719         computeLocalOffset(edgesToSend.size(), localOffset, nGlobalEdges);
00720 
00721         //Communicate the globalIDs, computed locally, to the other processors.
00722         for(ID index=0; index<edgesToSend.size(); index++)
00723             fEdgeToEdgeID[edgesToSend[index] ] = index + localOffset;
00724 
00725         allToAll(fEdgeToEdgeID, sendEdgesList_F, recvEdgesList_F);
00726 
00727 
00728         nEdges = edgeToFEdge.size();
00729         indexToEdgeID.resize(nEdges);
00730         for(int index=0; index<nEdges; index++)
00731            indexToEdgeID[index] = fEdgeToEdgeID[edgeToFEdge[index] ];
00732            
00733         //Compute vertices:
00734         std::vector<int> fCellsToSend; fCellsToSend.reserve(nCellsSolve_F);
00735 
00736 
00737         vertexToFCell.clear();
00738         vertexToFCell.reserve(nCells_F);
00739 
00740         fCellToVertex.assign(nCells_F,NotAnId);
00741         std::vector<int> fCellToVertexID(nCells_F, NotAnId);
00742 
00743         fCellsToReceive.clear();
00744 
00745        // if(! isDomainEmpty)
00746        // {
00747     fCellsToReceive.reserve(nCells_F-nCellsSolve_F);
00748     for(int i=0; i<nCells_F; i++)
00749     {
00750       bool isMine = i<nCellsSolve_F;
00751       bool belongsToLocalTriangle = false;
00752       bool belongsToAnyTriangle = false;
00753       int nEdg = nEdgesOnCells_F[i];
00754       for(int j=0; j<nEdg; j++)
00755       {
00756         ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F*i+j]-1);
00757         ID iTria = fVertexToTriangle[fVertex];
00758         ID triaId = fVertexToTriangleID[fVertex];
00759         belongsToLocalTriangle = belongsToLocalTriangle || (iTria != NotAnId);
00760         belongsToAnyTriangle = belongsToAnyTriangle || (triaId != NotAnId);
00761       }
00762 
00763       if( belongsToAnyTriangle && isMine)
00764       {
00765         fCellsToSend.push_back(i);
00766         if(!belongsToLocalTriangle)
00767           fCellsToReceive.push_back(i);
00768       }
00769 
00770       if( belongsToLocalTriangle)
00771       {
00772         fCellToVertex[i] = vertexToFCell.size();
00773         vertexToFCell.push_back(i);
00774       }
00775     }
00776 //        }
00777 
00778         //Compute the global number of vertices, and the localOffset on the local processor, such that a globalID = localOffset + index
00779         computeLocalOffset(fCellsToSend.size(), localOffset, nGlobalVertices);
00780 
00781         //Communicate the globalIDs, computed locally, to the other processors.
00782         for(int index=0; index< int(fCellsToSend.size()); index++)
00783             fCellToVertexID[fCellsToSend[index] ] = index + localOffset;
00784 
00785         allToAll(fCellToVertexID, sendCellsList_F, recvCellsList_F);
00786 
00787         nVertices = vertexToFCell.size();
00788         std::cout << "\n nvertices: " << nVertices << " " << nGlobalVertices << "\n"<<std::endl;
00789         indexToVertexID.resize(nVertices);
00790         for(int index=0; index<nVertices; index++)
00791             indexToVertexID[index] = fCellToVertexID[vertexToFCell[index] ];
00792         
00793         createReverseCellsExchangeLists(sendCellsListReversed, recvCellsListReversed, fVertexToTriangleID, fCellToVertexID);
00794 
00795 
00796         //construct the local vector vertices on triangles making sure the area is positive
00797         verticesOnTria.resize(nTriangles*3);
00798         double x[3], y[3], z[3];
00799         for(int index=0; index<nTriangles; index++)
00800         {
00801             int iTria = triangleToFVertex[index];
00802 
00803             for (int j=0; j<3; j++)
00804             {
00805               int iCell = cellsOnVertex_F[3*iTria+j]-1;
00806                 verticesOnTria[3*index+j] = fCellToVertex[iCell];
00807                 x[j] = xCell_F[iCell];
00808                 y[j] = yCell_F[iCell];
00809               //  z[j] = zCell_F[iCell];
00810             }
00811             if(signedTriangleArea(x,y) < 0)
00812               std::swap(verticesOnTria[3*index+1],verticesOnTria[3*index+2]);
00813         }
00814 
00815         //construct the local vector vertices on edges
00816         trianglesPositionsOnEdge.resize(2*nEdges);
00817         isVertexBoundary.assign(nVertices,false);
00818        
00819         verticesOnEdge.resize(2*nEdges);
00820 
00821 
00822         //contains the local id of a triangle and the global id of the edges of the triangle.
00823         //dataForGhostTria[4*i] contains the triangle id
00824         //dataForGhostTria[4*i+1+k] contains the global id of the edge (at position k = 0,1,2) of the triangle.
00825         //Possible Optimization: for our purposes it would be enough to store two of the three edges of a triangle.
00826         std::vector<int> dataForGhostTria(nVertices_F*4, NotAnId);
00827 
00828 
00829 //*
00830         for(int iV=0; iV<nVertices; iV++)
00831         {
00832       int fCell = vertexToFCell[iV];
00833       int nEdg = nEdgesOnCells_F[fCell];
00834       int j=0;
00835       bool isBoundary;
00836       do
00837       {
00838         int fVertex = verticesOnCell_F[maxNEdgesOnCell_F*fCell+j++]-1;
00839         isBoundary = !(verticesMask_F[fVertex] & 0x02);
00840       } while ((j<nEdg)&&(!isBoundary));
00841       isVertexBoundary[ iV  ] = isBoundary;
00842         }
00843  /*/
00844        for(int index=0; index<numBoundaryEdges; index++)
00845      {
00846       int fEdge = edgeToFEdge[index];
00847        int v1 = fCellToVertex[cellsOnEdge_F[2*fEdge]-1];
00848        int v2 = fCellToVertex[cellsOnEdge_F[2*fEdge+1]-1];
00849        isVertexBoundary[ v1 ] = isVertexBoundary[ v2 ] = true;
00850      }
00851 //*/
00852         //computing the position and local ids of the triangles adjacent to edges.
00853         //in the case an adjacent triangle is not local, the position and id will be NotAnId.
00854         //We will get the needed information about the non local edge leter on,
00855         //for this purpose we fill the vector dataForGhostTria.
00856         for(int index=0; index<nEdges; index++)
00857         {
00858             int p1,p2, v1,v2,v3, t, position;
00859             int fEdge = edgeToFEdge[index];
00860             int fCell1 = cellsOnEdge_F[2*fEdge]-1;
00861             int fCell2 = cellsOnEdge_F[2*fEdge+1]-1;
00862             verticesOnEdge[2*index] = p1 =fCellToVertex[fCell1];
00863             verticesOnEdge[2*index+1] = p2 =fCellToVertex[fCell2];
00864 
00865             for(int k=0;k<2 && (t = trianglesOnEdge[2*index+k]) != NotAnId; k++)
00866             {
00867                 v1 = verticesOnTria[3*t];
00868                 v2 = verticesOnTria[3*t+1];
00869                 v3 = verticesOnTria[3*t+2];
00870                 position = (((p1==v2)&&(p2==v3)) || ((p1==v3)&&(p2==v2))) + 2*(((p1==v1)&&(p2==v3)) || ((p1==v3)&&(p2==v1)));
00871                 trianglesPositionsOnEdge[2*index+k] = position;
00872                 int dataIndex = 4*triangleToFVertex[t];
00873 
00874 
00875                 dataForGhostTria[dataIndex] = t;
00876                 dataForGhostTria[dataIndex+position+1] = indexToEdgeID[index];
00877             }
00878         }
00879 
00880 
00881 
00882         createReverseEdgesExchangeLists(sendEdgesListReversed, recvEdgesListReversed, fVertexToTriangleID, fEdgeToEdgeID);
00883 
00884         //send the information about ghost elements.
00885         allToAll(dataForGhostTria, sendVerticesList_F, recvVerticesList_F, 4);
00886 
00887 
00888 
00889         //retrieving the position, local id and owner processor of non local triangles adjacent to interface edges.
00890         for(int index=0; index<(int)procOnInterfaceEdge.size(); index++)
00891         {
00892             int iEdge = procOnInterfaceEdge[index].first;
00893             int fEdge = edgeToFEdge[iEdge];
00894             ID fVertex1(verticesOnEdge_F[2*fEdge]-1), fVertex2(verticesOnEdge_F[2*fEdge+1]-1);
00895             if((ID)fVertexToTriangle[fVertex1] == NotAnId) fVertex2 = fVertex1;
00896             int dataIndex = 4*fVertex2;
00897             int edgeID = indexToEdgeID[iEdge];
00898             int position = (dataForGhostTria[dataIndex+2] == edgeID) + 2*(dataForGhostTria[dataIndex+3] == edgeID);
00899             trianglesOnEdge[2*iEdge+1] = dataForGhostTria[dataIndex];
00900             trianglesPositionsOnEdge[2*iEdge+1] = position;
00901          }
00902 
00903         if( isDomainEmpty )
00904             return;
00905 
00906    //     ArrayRCP<RCP<Albany::MeshSpecsStruct> > meshSpecs =  discFactory.createMeshSpecs();
00907 
00908         //construct the local vector of coordinates
00909   //      std::vector<double> verticesCoords(3*nVertices);
00910 
00911 
00912  //   for(int index=0; index<nVertices; index++)
00913  //   {
00914  //     int iCell = vertexToFCell[index];
00915  //     verticesCoords[index*3] = xCell_F[iCell]/unit_length;
00916  //     verticesCoords[index*3 + 1] = yCell_F[iCell]/unit_length;
00917  //     verticesCoords[index*3 + 2] = zCell_F[iCell]/unit_length;
00918  //   }
00919 
00920     mpiComm = Albany::createEpetraCommFromMpiComm(reducedComm);
00921 //                std::string xmlfilename = "albany_input.xml";
00922 
00923 // GET slvrfctry STUFF FROM 3D GRID BELOW
00924 
00925 //                meshStruct2D = Teuchos::rcp(new Albany::MpasSTKMeshStruct(discParams, mpiComm, indexToTriangleID, verticesOnTria, nGlobalTriangles));
00926 //                meshStruct2D->constructMesh(mpiComm, discParams, sis, indexToVertexID, verticesCoords, isVertexBoundary, nGlobalVertices,
00927 //                           verticesOnTria, isBoundaryEdge, trianglesOnEdge, trianglesPositionsOnEdge,
00928 //                           verticesOnEdge, indexToEdgeID, nGlobalEdges, meshStruct->getMeshSpecs()[0]->worksetSize);
00929 
00930         /*
00931         //initialize the mesh
00932         iceProblemPtr->mesh2DPtr.reset(new RegionMesh<LinearTriangle>() );
00933 
00934         //construct the mesh nodes
00935         constructNodes( *(iceProblemPtr->mesh2DPtr), indexToVertexID, verticesCoords, isVertexBoundary, nGlobalVertices, 3);
00936         
00937         //construct the mesh elements
00938         constructElements( *(iceProblemPtr->mesh2DPtr), indexToTriangleID, verticesOnTria, nGlobalTriangles);
00939 
00940         //construct the mesh facets
00941         constructFacets( *(iceProblemPtr->mesh2DPtr), isBoundaryEdge, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge, indexToEdgeID, procOnInterfaceEdge, nGlobalEdges, 3);
00942 
00943         Switch sw;
00944         std::vector<bool> elSign;
00945         checkVolumes( *(iceProblemPtr->mesh2DPtr),elSign, sw );
00946 
00947         */
00948   }
00949 
00950 
00951 
00952     void velocity_solver_extrude_3d_grid(double const * levelsRatio_F, double const * lowerSurface_F, double const * thickness_F)
00953     {
00954         if(isDomainEmpty)
00955             return;
00956 
00957         layersRatio.resize(nLayers);
00958         // !!Indexing of layers is reversed
00959         for(int i=0; i<nLayers; i++)
00960             layersRatio[i] = levelsRatio_F[nLayers-1-i];
00961         //std::copy(levelsRatio_F, levelsRatio_F+nLayers, layersRatio.begin());
00962 
00963         levelsNormalizedThickness.resize(nLayers+1);
00964 
00965         levelsNormalizedThickness[0] = 0;
00966         for(int i=0; i< nLayers; i++)
00967             levelsNormalizedThickness[i+1] = levelsNormalizedThickness[i]+layersRatio[i];
00968 
00969     std::vector<int> mpasIndexToVertexID(nVertices);
00970     for(int i=0; i<nVertices; i++)
00971       mpasIndexToVertexID[i] = indexToCellID_F[vertexToFCell[i]];
00972 
00973 
00974     //construct the local vector of coordinates
00975       std::vector<double> verticesCoords(3*nVertices);
00976 
00977 
00978       for(int index=0; index<nVertices; index++)
00979       {
00980         int iCell = vertexToFCell[index];
00981         verticesCoords[index*3] = xCell_F[iCell]/unit_length;
00982         verticesCoords[index*3 + 1] = yCell_F[iCell]/unit_length;
00983         verticesCoords[index*3 + 2] = zCell_F[iCell]/unit_length;
00984       }
00985 
00986 
00987       slvrfctry = Teuchos::rcp(new Albany::SolverFactory("albany_input.xml", reducedComm));
00988       discParams = Teuchos::sublist(Teuchos::rcp(&slvrfctry->getParameters(),false), "Discretization", true);
00989       Teuchos::RCP<Albany::StateInfoStruct> sis=Teuchos::rcp(new Albany::StateInfoStruct);
00990 /*
00991       meshStruct = Teuchos::rcp(new Albany::MpasSTKMeshStruct(discParams, mpiComm, indexToTriangleID, verticesOnTria, nGlobalTriangles,nLayers,Ordering));
00992       meshStruct->constructMesh(mpiComm, discParams, sis, indexToVertexID, verticesCoords, isVertexBoundary, nGlobalVertices,
00993                  verticesOnTria, isBoundaryEdge, trianglesOnEdge, trianglesPositionsOnEdge,
00994                  verticesOnEdge, indexToEdgeID, nGlobalEdges, indexToTriangleID, meshStruct->getMeshSpecs()[0]->worksetSize,nLayers,Ordering);
00995   */
00996       Albany::AbstractFieldContainer::FieldContainerRequirements req;
00997       req.push_back("Surface Height");
00998       req.push_back("Temperature");
00999       req.push_back("Basal Friction");
01000       req.push_back("Thickness");
01001       int neq=2;
01002       meshStruct = Teuchos::rcp(new Albany::MpasSTKMeshStruct(discParams, mpiComm, indexToTriangleID, nGlobalTriangles,nLayers,Ordering));
01003       meshStruct->constructMesh(mpiComm, discParams, neq, req, sis, indexToVertexID, mpasIndexToVertexID, verticesCoords, isVertexBoundary, nGlobalVertices,
01004                 verticesOnTria, isBoundaryEdge, trianglesOnEdge, trianglesPositionsOnEdge,
01005                 verticesOnEdge, indexToEdgeID, nGlobalEdges, indexToTriangleID, meshStruct->getMeshSpecs()[0]->worksetSize,nLayers,Ordering);
01006     }
01007 }
01008 
01009 
01010 
01011   void get_tetraP1_velocity_on_FEdges(double * uNormal, const std::vector<double>& velocityOnVertices, const std::vector<int>& edgeToFEdge, const std::vector<int>& mpasIndexToVertexID)
01012   {
01013 
01014     int columnShift = (Ordering == 1) ? 1 : nVertices;
01015     int layerShift = (Ordering == 0) ? 1 : nLayers + 1;
01016 
01017     UInt nPoints3D = nVertices * (nLayers + 1);
01018 
01019     //the velocity on boundary edges is set to zero by construction
01020     for(int i=numBoundaryEdges; i< nEdges; i++)
01021     {
01022        ID lId0 = verticesOnEdge[2*i];
01023        ID lId1 = verticesOnEdge[2*i+1];
01024        int iCell0 = vertexToFCell[lId0];
01025        int iCell1 = vertexToFCell[lId1];
01026 
01027        double nx = xCell_F[iCell1] - xCell_F[iCell0];
01028        double ny = yCell_F[iCell1] - yCell_F[iCell0];
01029        double n = sqrt(nx*nx+ny*ny);
01030        nx /= n;
01031        ny /= n;
01032 
01033        int iEdge = edgeToFEdge[i];
01034        //prism lateral face is splitted with the diagonal that start from p0 (by construction).
01035        double coeff0 = (mpasIndexToVertexID[lId0] > mpasIndexToVertexID[lId1]) ?  1./3. : 1./6.;
01036        double coeff1 = 0.5 - coeff0;
01037        for(int il=0; il < nLayers; il++)
01038        {
01039          int ilReversed = nLayers-il-1;
01040          ID lId3D1 = il * columnShift + layerShift * lId0;
01041          ID lId3D2 = il * columnShift + layerShift * lId1;
01042          //not accurate
01043          uNormal[iEdge*nLayers+ilReversed] = nx*(coeff1*velocityOnVertices[lId3D1] + coeff0*velocityOnVertices[lId3D2] + coeff0*velocityOnVertices[lId3D1+columnShift] + coeff1*velocityOnVertices[lId3D2+columnShift]);
01044 
01045          lId3D1 += nPoints3D;
01046          lId3D2 += nPoints3D;
01047          uNormal[iEdge*nLayers+ilReversed] += ny*(coeff1*velocityOnVertices[lId3D1] + coeff0*velocityOnVertices[lId3D2] + coeff0*velocityOnVertices[lId3D1+columnShift] + coeff1*velocityOnVertices[lId3D2+columnShift]);
01048        }
01049     }
01050 
01051   }
01052 
01053   void get_prism_velocity_on_FEdges(double * uNormal, const std::vector<double>& velocityOnVertices, const std::vector<int>& edgeToFEdge)
01054   {
01055 
01056       //fix this
01057     int columnShift = (Ordering == 1) ? 1 : nVertices;
01058     int layerShift = (Ordering == 0) ? 1 : nLayers + 1;
01059 
01060     UInt nPoints3D = nVertices * (nLayers + 1);
01061 
01062     for(int i=numBoundaryEdges; i< nEdges; i++)
01063     {
01064       ID lId0 = verticesOnEdge[2*i];
01065       ID lId1 = verticesOnEdge[2*i+1];
01066       int iCell0 = vertexToFCell[lId0];
01067       int iCell1 = vertexToFCell[lId1];
01068 
01069       double nx = xCell_F[iCell1] - xCell_F[iCell0];
01070       double ny = yCell_F[iCell1] - yCell_F[iCell0];
01071       double n = sqrt(nx*nx+ny*ny);
01072       nx /= n;
01073       ny /= n;
01074 
01075        ID iEdge = edgeToFEdge[i];
01076 
01077        for(int il=0; il < nLayers; il++)
01078        {
01079          int ilReversed = nLayers-il-1;
01080          ID lId3D1 = il * columnShift + layerShift * lId0;
01081          ID lId3D2 = il * columnShift + layerShift * lId1;
01082          //not accurate
01083          uNormal[iEdge*nLayers+ilReversed] = 0.25*nx*(velocityOnVertices[lId3D1] + velocityOnVertices[lId3D2] + velocityOnVertices[lId3D1+columnShift] + velocityOnVertices[lId3D2+columnShift]);
01084 
01085          lId3D1 += nPoints3D;
01086          lId3D2 += nPoints3D;
01087          uNormal[iEdge*nLayers+ilReversed] += 0.25*ny*(velocityOnVertices[lId3D1] + velocityOnVertices[lId3D2] + velocityOnVertices[lId3D1+columnShift] + velocityOnVertices[lId3D2+columnShift]);
01088        }
01089     }
01090 
01091   }
01092 
01093 
01094     void mapVerticesToCells(const std::vector<double>& velocityOnVertices, double* velocityOnCells, int fieldDim, int numLayers, int ordering)
01095     {
01096       int lVertexColumnShift = (ordering == 1) ? 1 : nVertices;
01097     int vertexLayerShift = (ordering == 0) ? 1 : numLayers+1;
01098 
01099     int nVertices3D = nVertices*(numLayers+1);
01100     for ( UInt j = 0 ; j < nVertices3D ; ++j )
01101     {
01102        int ib = (ordering == 0)*(j%lVertexColumnShift) + (ordering == 1)*(j/vertexLayerShift);
01103        int il = (ordering == 0)*(j/lVertexColumnShift) + (ordering == 1)*(j%vertexLayerShift);
01104 
01105        int iCell = vertexToFCell[ ib ];
01106        int cellIndex = iCell * (numLayers+1) + il;
01107        int vertexIndex = j;
01108        for( int dim = 0; dim < fieldDim; dim++ )
01109        {
01110          velocityOnCells[cellIndex] = velocityOnVertices[vertexIndex];
01111          cellIndex += nCells_F * (numLayers+1);
01112          vertexIndex += nVertices3D;
01113        }
01114     }
01115 
01116         for( int dim = 0; dim < fieldDim; dim++ )
01117         {
01118           allToAll(&velocityOnCells[dim*nCells_F*(numLayers+1)],  &sendCellsListReversed,&recvCellsListReversed, (numLayers+1));
01119           allToAll(&velocityOnCells[dim*nCells_F*(numLayers+1)],  sendCellsList_F, recvCellsList_F, (numLayers+1));
01120         }
01121     }
01122 
01123 
01124     void createReverseCellsExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, const std::vector<int>& fVertexToTriangleID, const std::vector<int>& fCellToVertexID)
01125     {
01126         sendListReverse_F.clear();
01127         receiveListReverse_F.clear();
01128         //std::map<int, std::vector<int> > sendMap;
01129         std::map<int, std::map<int, int> > sendMap, receiveMap;
01130         std::vector<int> cellsProcId(nCells_F), trianglesProcIds(nVertices_F);
01131         getProcIds(cellsProcId, recvCellsList_F);
01132         getProcIds(trianglesProcIds, recvVerticesList_F);
01133 
01134         //std::cout << "SendList " ;
01135         for(int i=0; i< nVertices; i++)
01136         {
01137             int iCell = vertexToFCell[i];
01138             if(iCell<nCellsSolve_F) continue;
01139             bool belongToTriaOnSameProc = false;
01140             int j(0);
01141             int nEdg = nEdgesOnCells_F[iCell];
01142             do
01143             {
01144                 ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F*iCell+j]-1);
01145                 ID triaId = fVertexToTriangleID[fVertex];
01146                 belongToTriaOnSameProc = (triaId != NotAnId)&&(trianglesProcIds[fVertex] == cellsProcId[iCell]);
01147             }while( (belongToTriaOnSameProc==false) && (++j < nEdg));
01148             if(!belongToTriaOnSameProc)
01149             {
01150             sendMap[cellsProcId[iCell]].insert(std::make_pair(fCellToVertexID[iCell], iCell));
01151            // std::cout<< "(" << cellsProcId[iCell] << "," << iCell << ") ";
01152             }
01153 
01154         }
01155         //std::cout <<std::endl;
01156 
01157         for(std::map<int, std::map<int,int>  >::const_iterator it=sendMap.begin(); it!=sendMap.end(); it++)
01158     {
01159       std::vector<int> sendVec(it->second.size());
01160       int i=0;
01161       for(std::map<int,int>::const_iterator iter= it->second.begin(); iter != it->second.end(); iter++)
01162         sendVec[i++] = iter->second;
01163       sendListReverse_F.push_back(exchange(it->first, &sendVec[0], &sendVec[0]+sendVec.size()));
01164     }
01165 
01166         //std::cout << "ReceiveList " ;
01167         for(UInt i=0; i< fCellsToReceive.size(); i++)
01168         {
01169             int iCell = fCellsToReceive[i];
01170             int nEdg = nEdgesOnCells_F[iCell];
01171             for(int j=0; j<nEdg; j++)
01172             {
01173                 ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F*iCell+j]-1);
01174                 ID triaId = fVertexToTriangleID[fVertex];
01175                 if(triaId != NotAnId)
01176                 {
01177                     receiveMap[trianglesProcIds[fVertex] ].insert(std::make_pair(fCellToVertexID[iCell], iCell));
01178                    // std::cout<< "(" << trianglesProcIds[fVertex] << "," << iCell << ") ";
01179                 }
01180             }
01181         }
01182         //std::cout <<std::endl;
01183 
01184         for(std::map<int, std::map<int,int> >::const_iterator it=receiveMap.begin(); it!=receiveMap.end(); it++)
01185     {
01186       std::vector<int> receiveVec(it->second.size());
01187       int i=0;
01188       for(std::map<int,int>::const_iterator iter= it->second.begin(); iter != it->second.end(); iter++)
01189         receiveVec[i++] = iter->second;
01190       receiveListReverse_F.push_back(exchange(it->first, &receiveVec[0], &receiveVec[0]+receiveVec.size()));
01191     }
01192     }
01193 
01194     void createReverseEdgesExchangeLists(exchangeList_Type& sendListReverse_F, exchangeList_Type& receiveListReverse_F, const std::vector<int>& fVertexToTriangleID, const std::vector<int>& fEdgeToEdgeID)
01195     {
01196         sendListReverse_F.clear();
01197         receiveListReverse_F.clear();
01198         //std::map<int, std::vector<int> > sendMap;
01199         std::map<int, std::map<int, int> > sendMap, receiveMap;
01200         std::vector<int> edgesProcId(nEdges_F), trianglesProcIds(nVertices_F);
01201         getProcIds(edgesProcId, recvEdgesList_F);
01202         getProcIds(trianglesProcIds, recvVerticesList_F);
01203 
01204         //std::cout << "EdgesSendList " ;
01205         for(int i=0; i< nEdges; i++)
01206         {
01207             int iEdge = edgeToFEdge[i];
01208             if(iEdge<nEdgesSolve_F) continue;
01209             bool belongToTriaOnSameProc = false;
01210             int j(0);
01211             do
01212             {
01213                 ID fVertex(verticesOnEdge_F[2*iEdge+j]-1);
01214                 ID triaId = fVertexToTriangleID[fVertex];
01215                 belongToTriaOnSameProc = (triaId != NotAnId)&&(trianglesProcIds[fVertex] == edgesProcId[iEdge]);
01216             }while( (belongToTriaOnSameProc==false) && (++j < 2));
01217             if(!belongToTriaOnSameProc){
01218                 sendMap[edgesProcId[iEdge]].insert(std::make_pair( fEdgeToEdgeID[iEdge], iEdge));
01219                 //std::cout<< "(" << edgesProcId[iEdge] << "," << iEdge << ") ";
01220             }
01221         }
01222         //std::cout <<std::endl;
01223 
01224         for(std::map<int, std::map<int,int>  >::const_iterator it=sendMap.begin(); it!=sendMap.end(); it++)
01225         {
01226           std::vector<int> sendVec(it->second.size());
01227           int i=0;
01228           for(std::map<int,int>::const_iterator iter= it->second.begin(); iter != it->second.end(); iter++)
01229             sendVec[i++] = iter->second;
01230             sendListReverse_F.push_back(exchange(it->first, &sendVec[0], &sendVec[0]+sendVec.size()));
01231         }
01232 
01233         //std::cout << "EdgesReceiveList " ;
01234         for(UInt i=0; i< edgesToReceive.size(); i++)
01235         {
01236             int iEdge = edgesToReceive[i];
01237             for(int j=0; j<2; j++)
01238             {
01239                 ID fVertex(verticesOnEdge_F[2*iEdge+j]-1);
01240                 ID triaId = fVertexToTriangleID[fVertex];
01241                 if(triaId != NotAnId){
01242                     receiveMap[trianglesProcIds[fVertex] ].insert(std::make_pair( fEdgeToEdgeID[iEdge], iEdge));
01243                   //  std::cout<< "(" << trianglesProcIds[fVertex] << "," << iEdge << ") ";
01244                 }
01245             }
01246         }
01247        // std::cout <<std::endl;
01248 
01249         for(std::map<int, std::map<int,int> >::const_iterator it=receiveMap.begin(); it!=receiveMap.end(); it++)
01250         {
01251           std::vector<int> receiveVec(it->second.size());
01252           int i=0;
01253           for(std::map<int,int>::const_iterator iter= it->second.begin(); iter != it->second.end(); iter++)
01254             receiveVec[i++] = iter->second;
01255             receiveListReverse_F.push_back(exchange(it->first, &receiveVec[0], &receiveVec[0]+receiveVec.size()));
01256         }
01257     }
01258 
01259 /*
01260     void testAllToAll(const std::vector<int>& fCellToVertex, const std::vector<int>& fCellToVertexID, const std::vector<int>& fVertexToTriangleID )
01261     {
01262 
01263      //   exchangeList_Type sendList, receiveList;
01264      //   createReverseCellsExchangeLists(sendList, receiveList, fVertexToTriangleID);
01265         std::vector<int> vertexIDOnCells(nCells_F, NotAnId), verticesIDs(nVertices, NotAnId);
01266 
01267 /*
01268         std::map<int, std::vector<int> > sendMap, receiveMap;
01269         std::vector<int> cellsProcId(nCells_F), trianglesProcIds(nVertices_F);
01270         getProcIds(cellsProcId, recvCellsList_F);
01271         getProcIds(trianglesProcIds, recvVerticesList_F);
01272 
01273         for(UInt i=0; i< nVertices; i++)
01274         {
01275             int iCell = vertexToFCell[i];
01276             if(iCell<nCellsSolve_F) continue;
01277             sendMap[cellsProcId[iCell]].push_back(iCell);
01278         }
01279 
01280         for(std::map<int, std::vector<int> >::const_iterator it=sendMap.begin(); it!=sendMap.end(); it++)
01281             sendList.push_back(exchange(it->first, &it->second[0], &it->second[0]+it->second.size()));
01282 
01283         for(UInt i=0; i< fCellsToReceive.size(); i++)
01284         {
01285             int iCell = fCellsToReceive[i];
01286             int nEdg = nEdgesOnCells_F[iCell];
01287             for(int j=0; j<nEdg; j++)
01288             {
01289                 ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F*iCell+j]-1);
01290                 ID triaId = fVertexToTriangleID[fVertex];
01291                 if(triaId != NotAnId)
01292                     receiveMap[trianglesProcIds[fVertex] ].push_back(iCell);
01293             }
01294 
01295         }
01296 
01297         for(std::map<int, std::vector<int> >::const_iterator it=receiveMap.begin(); it!=receiveMap.end(); it++)
01298             receiveList.push_back(exchange(it->first, &it->second[0], &it->second[0]+it->second.size()));
01299 
01300 
01301 
01302 
01303 
01304 
01305 
01306         int me;
01307         MPI_Comm_rank ( comm, &me );
01308         if(me == 8)
01309             vertexIDOnCells[256] = 1;
01310 
01311         allToAll(vertexIDOnCells,  sendCellsList_F, recvCellsList_F, 1);
01312 
01313         if(me == 8)
01314             vertexIDOnCells[256] = NotAnId;
01315 
01316  //       for(int index=0; index<nVertices; index++)
01317  //       {
01318  //           if( vertexIDOnCells[vertexToFCell[ index ]] == 1)
01319  //              std::cout<< "Arrivato! "<<vertexToFCell[ index ]  <<std::endl;
01320  //       }
01321 
01322         for(int index=0; index<nCells_F; index++)
01323                 {
01324                     if( vertexIDOnCells[ index ] == 1)
01325                         std::cout<< "Arrivato! "<<index   <<std::endl;
01326                 }
01327 
01328 
01329 
01330 
01331         allToAll(vertexIDOnCells,  &sendList, &receiveList,1);
01332 
01333 
01334         if((me == 8)&&(vertexIDOnCells[256] ==1))
01335                 std::cout<< "Ritornato!"<<std::endl;
01336 
01337         if(me == 8)
01338                     vertexIDOnCells[256] = NotAnId;
01339 
01340 
01341 
01342 
01343 * /
01344 
01345 
01346         for(int index=0; index<nVertices; index++)
01347             vertexIDOnCells[vertexToFCell[ index ]] = indexToVertexID[index];
01348 
01349      //   std::vector<int> tmp_vec(vertexIDOnCells);
01350 
01351 
01352         allToAll(vertexIDOnCells, &sendCellsListReversed, &recvCellsListReversed,1);
01353 
01354 //        for(UInt i=0; i< fCellsToReceive.size(); i++)
01355  //       {
01356   //          UInt iCell = fCellsToReceive[i];
01357    //         std::cout<< "received " << fCellToVertexID[iCell] << " vertex: " <<  tmp_vec[iCell] <<", it was " << vertexIDOnCells[iCell] << ", cells: " << fCellToVertex[iCell] << "<" <<nVertices << std::endl;
01358          //   vertexIDOnCells[iCell] = tmp_vec[iCell];
01359    //      }
01360 
01361         allToAll(vertexIDOnCells,  sendCellsList_F, recvCellsList_F, 1);
01362 
01363         for(int index=0; index<nVertices; index++)
01364             verticesIDs[index] = vertexIDOnCells[vertexToFCell[ index ]];
01365 
01366 
01367 
01368         double sum(0), total(0);
01369         for(UInt i=0; i<nVertices; i++)
01370         {
01371             if(verticesIDs[i]!=indexToVertexID[i])
01372             {
01373                 std::cout << "ecco: " << indexToVertexID[i] << "!=" << verticesIDs[i] <<  ", isMine? " << (vertexToFCell[i] < nCellsSolve_F) <<std::endl;
01374             }
01375         //    sum += fabs(pippo[i]-velocityOnVertices[i]);
01376         }
01377            // sum += fabs(pippo[i]-velocityOnVertices[i]);
01378 
01379       //  MPI_Reduce ( &sum, &total, 1, MPI_DOUBLE, MPI_SUM, 0, reducedComm );
01380 
01381 //        std::cout << "Error: " << sum << " " << total << std::endl;
01382 
01383      //   exit(0);
01384 
01385 
01386     }
01387 */
01388     void mapCellsToVertices(const std::vector<double>& velocityOnCells, std::vector<double>& velocityOnVertices, int fieldDim, int numLayers, int ordering)
01389     {
01390       int lVertexColumnShift = (ordering == 1) ? 1 : nVertices;
01391     int vertexLayerShift = (ordering == 0) ? 1 : numLayers+1;
01392 
01393     int nVertices3D = nVertices*(numLayers+1);
01394     for ( UInt j = 0 ; j < nVertices3D ; ++j )
01395     {
01396        int ib = (ordering == 0)*(j%lVertexColumnShift) + (ordering == 1)*(j/vertexLayerShift);
01397        int il = (ordering == 0)*(j/lVertexColumnShift) + (ordering == 1)*(j%vertexLayerShift);
01398 
01399        int iCell = vertexToFCell[ ib ];
01400        int cellIndex = iCell * (numLayers+1) + il;
01401        int vertexIndex = j;
01402        for( int dim = 0; dim < fieldDim; dim++ )
01403        {
01404          velocityOnVertices[vertexIndex] = velocityOnCells[cellIndex];
01405          cellIndex += nCells_F * (numLayers+1);
01406          vertexIndex += nVertices3D;
01407        }
01408     }
01409     }
01410 
01411     bool isGhostTriangle(int i, double relTol)
01412     {
01413       double x[3], y[3], area;
01414 
01415     for(int j=0; j<3; j++)
01416     {
01417         int iCell = cellsOnVertex_F[3*i+j]-1;
01418         x[j] = xCell_F[iCell];
01419         y[j] = yCell_F[iCell];
01420     }
01421 
01422     area = std::fabs(signedTriangleArea(x,y));
01423     return false; //(std::fabs(areaTriangle_F[i]-area)/areaTriangle_F[i] > relTol);
01424     }
01425 
01426     double signedTriangleArea(const double* x, const double* y)
01427     {
01428       double u[2] = {x[1]-x[0], y[1]-y[0]};
01429       double v[2] = {x[2]-x[0], y[2]-y[0]};
01430 
01431       return 0.5*(u[0]*v[1]-u[1]*v[0]);
01432     }
01433 
01434     double signedTriangleAreaOnSphere(const double* x, const double* y, const double *z)
01435   {
01436     double u[3] = {x[1]-x[0], y[1]-y[0], z[1]-z[0]};
01437     double v[3] = {x[2]-x[0], y[2]-y[0], z[2]-z[0]};
01438 
01439     double crossProduct[3] = {u[1]*v[2]-u[2]*v[1], u[2]*v[0]-u[0]*v[2], u[0]*v[1]-u[1]*v[0]};
01440     double area = 0.5*std::sqrt(crossProduct[0]*crossProduct[0]+crossProduct[1]*crossProduct[1]+crossProduct[2]*crossProduct[2]);
01441     return (crossProduct[0]*x[0]+crossProduct[1]*y[0]+crossProduct[2]*z[0] > 0) ? area : -area;
01442   }
01443 
01444 
01445     //TO BE FIXED, Access To verticesOnCell_F is not correct
01446     void extendMaskByOneLayer(int const* verticesMask_F, std::vector<int>& extendedFVerticesMask)
01447     {
01448       extendedFVerticesMask.resize(nVertices_F);
01449       extendedFVerticesMask.assign(&verticesMask_F[0], &verticesMask_F[0] + nVertices_F);
01450       for( int i=0; i<nCells_F; i++)
01451     {
01452         bool belongsToMarkedTriangle=false;
01453         int nEdg = nEdgesOnCells_F[i];
01454         for(UInt k=0; k<nEdg && !belongsToMarkedTriangle; k++)
01455           belongsToMarkedTriangle = belongsToMarkedTriangle || verticesMask_F[verticesOnCell_F[maxNEdgesOnCell_F*i+k]-1];
01456         if(belongsToMarkedTriangle)
01457           for(UInt k=0; k<nEdg; k++)
01458           {
01459             ID fVertex(verticesOnCell_F[maxNEdgesOnCell_F*i+k]-1);
01460             extendedFVerticesMask[fVertex] = !isGhostTriangle(fVertex);
01461           }
01462     }
01463     }
01464 
01465 
01466     void import2DFields(double const * lowerSurface_F, double const * thickness_F,
01467                       double const * beta_F, double eps)
01468     {
01469       elevationData.assign(nVertices, 1e10);
01470     thicknessData.assign(nVertices, 1e10);
01471     std::map<int,int> bdExtensionMap;
01472     if(beta_F != 0)
01473       betaData.assign(nVertices,1e10);
01474 
01475     for(int iV=0; iV<nVertices; iV++)
01476     {
01477       if (isVertexBoundary[iV])
01478       {
01479         int c;
01480         int fCell = vertexToFCell[iV];
01481         int nEdg = nEdgesOnCells_F[fCell];
01482         for(int j=0; j<nEdg; j++)
01483         {
01484           int fEdge = edgesOnCell_F[maxNEdgesOnCell_F*fCell+j]-1;
01485           bool keep = (mask[verticesOnEdge_F[2*fEdge]-1]& 0x02) && (mask[verticesOnEdge_F[2*fEdge+1]-1]& 0x02);
01486                 if(!keep) continue;
01487 
01488                 int c0 = cellsOnEdge_F[2*fEdge]-1;
01489                 int c1 = cellsOnEdge_F[2*fEdge+1]-1;
01490                 c = (fCellToVertex[c0] == iV) ? c1 : c0;
01491                 double elev = thickness_F[c]+lowerSurface_F[c];// - 1e-8*std::sqrt(pow(xCell_F[c0],2)+std::pow(yCell_F[c0],2));
01492           if(elevationData[ iV ] > elev){
01493             elevationData[ iV ] = elev;
01494             bdExtensionMap[ iV ] = c;
01495           }
01496         }
01497       }
01498     }
01499 
01500     for(std::map<int,int>::iterator it=bdExtensionMap.begin(); it != bdExtensionMap.end(); ++it)
01501     {
01502       int iv = it->first;
01503       int ic = it->second;
01504       thicknessData[ iv ] = thickness_F[ic]/unit_length+eps;//- 1e-8*std::sqrt(pow(xCell_F[ic],2)+std::pow(yCell_F[ic],2));
01505       elevationData[ iv ] = std::max(thicknessData[ iv ] + lowerSurface_F[ic]/unit_length, 118./1028.*thicknessData[ iv ]);
01506       if(beta_F != 0)
01507       //  betaData[ iv ] = beta_F[ic]/unit_length;
01508           betaData[ iv ] = (lowerSurface_F[ic]> -910./1028.*thickness_F[ic]) ? beta_F[ic]/unit_length : 0;
01509     }
01510 
01511     for(int index=0; index<nVertices; index++)
01512     {
01513       int iCell = vertexToFCell[index];
01514 
01515       if (!isVertexBoundary[index])
01516       {
01517         thicknessData[ index ] = thickness_F[iCell]/unit_length+eps;
01518         elevationData[ index ] = std::max((lowerSurface_F[iCell]/unit_length)+thicknessData[ index ], 118./1028.*thicknessData[ index ]);
01519       }
01520     }
01521 
01522     if(beta_F != 0)
01523     {
01524       for(int index=0; index<nVertices; index++)
01525       {
01526         int iCell = vertexToFCell[index];
01527 
01528         if (!isVertexBoundary[index])
01529           //betaData[ index ] = beta_F[iCell]/unit_length;
01530             betaData[ index ] = (lowerSurface_F[iCell]> -910./1028.*thickness_F[iCell]) ? beta_F[iCell]/unit_length : 0;
01531       }
01532     }
01533 
01534     }
01535 
01536     void importP0Temperature(double const * temperature_F)
01537     {
01538       int lElemColumnShift = (Ordering == 1) ? 3 : 3*indexToTriangleID.size();
01539       int elemLayerShift = (Ordering == 0) ? 3 : 3*nLayers;
01540         temperatureOnTetra.resize(3*nLayers*indexToTriangleID.size());
01541         for(int index=0; index<nTriangles; index++)
01542         {
01543             for(int il=0; il < nLayers; il++)
01544             {
01545                 double temperature = 0;
01546                 int ilReversed = nLayers-il-1;
01547                 int nPoints=0;
01548                 for(int iVertex=0; iVertex<3; iVertex++)
01549                 {
01550                   int v = verticesOnTria[iVertex + 3*index];
01551                   if(!isVertexBoundary[v])
01552                   {
01553                     int iCell = vertexToFCell[ v ];
01554                     temperature  += temperature_F[iCell * nLayers + ilReversed];
01555                     nPoints++;
01556                   }
01557                 }
01558                 if(nPoints==0)
01559                   temperature = T0;
01560                 else
01561                   temperature = temperature/nPoints + T0;
01562                 for(int k=0; k<3; k++)
01563                     temperatureOnTetra[index * elemLayerShift + il*lElemColumnShift + k] = temperature;
01564             }
01565         }
01566 
01567     }
01568 
01569 
01570 
01571   void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id)
01572   {
01573     int numProcs, me;
01574     MPI_Group world_group_id, reduced_group_id;
01575     MPI_Comm_size ( comm, &numProcs );
01576     MPI_Comm_rank ( comm, &me );
01577     std::vector<int> haveElements(numProcs);
01578     int nonEmpty = int(nLocalEntities > 0);
01579     MPI_Allgather( &nonEmpty, 1, MPI_INT, &haveElements[0], 1, MPI_INT, comm);
01580     std::vector<int> ranks;
01581     for (int i=0; i<numProcs; i++)
01582     {
01583       if(haveElements[i])
01584         ranks.push_back(i);
01585     }
01586 
01587     MPI_Comm_group( comm, &world_group_id );
01588     MPI_Group_incl ( world_group_id, ranks.size(), &ranks[0], &reduced_group_id );
01589         MPI_Comm_create ( comm, reduced_group_id, &reduced_comm_id );
01590   }
01591 
01592 
01593   void computeLocalOffset(int nLocalEntities, int& localOffset, int& nGlobalEntities)
01594   {
01595     int numProcs, me;
01596     MPI_Comm_size ( comm, &numProcs );
01597     MPI_Comm_rank ( comm, &me );
01598     std::vector<int> offsetVec(numProcs);
01599 
01600     MPI_Allgather( &nLocalEntities, 1, MPI_INT, &offsetVec[0], 1, MPI_INT, comm);
01601 
01602     localOffset = 0;
01603     for(int i = 0; i<me; i++)
01604       localOffset += offsetVec[i];
01605 
01606     nGlobalEntities = localOffset;
01607     for(int i = me; i< numProcs; i++)
01608       nGlobalEntities += offsetVec[i];
01609   }
01610 
01611 
01612   void getProcIds(std::vector<int>& field,int const * recvArray)
01613   {
01614       int me;
01615       MPI_Comm_rank ( comm, &me );
01616       field.assign(field.size(),me);
01617 
01618         //unpack recvArray and set the proc rank into filed
01619         for(int i(1), procID, size; i< recvArray[0]; i += size)
01620         {
01621             procID = recvArray[i++];
01622             size = recvArray[i++];
01623             if(procID == me) continue;
01624             for(int k=i; k<i+size; k++)
01625                 field[recvArray[k]] = procID;
01626         }
01627   }
01628 
01629   void getProcIds(std::vector<int>& field, exchangeList_Type const * recvList)
01630     {
01631         int me;
01632         MPI_Comm_rank ( comm, &me );
01633         field.assign(field.size(),me);
01634         exchangeList_Type::const_iterator it;
01635 
01636         for(it = recvList->begin(); it != recvList->end(); ++it )
01637         {
01638             if(it->procID == me) continue;
01639             for(int k=0; k< (int)it->vec.size(); k++)
01640                   field[it->vec[k]] = it->procID;
01641         }
01642     }
01643 
01644 
01645   exchangeList_Type unpackMpiArray(int const * array)
01646   {
01647       exchangeList_Type list;
01648       for(int i(1), procID, size; i< array[0]; i += size)
01649         {
01650             procID = array[i++];
01651             size = array[i++];
01652             list.push_back(exchange(procID, &array[i], &array[i+size]) );
01653         }
01654       return list;
01655   }
01656 
01657 
01658   void allToAll(std::vector<int>& field, int const * sendArray, int const * recvArray, int fieldDim)
01659   {
01660       exchangeList_Type sendList, recvList;
01661 
01662     //unpack sendArray and build the sendList class
01663     for(int i(1), procID, size; i< sendArray[0]; i += size)
01664     {
01665       procID = sendArray[i++];
01666       size = sendArray[i++];
01667       sendList.push_back(exchange(procID, &sendArray[i], &sendArray[i+size], fieldDim) );
01668     }
01669 
01670     //unpack recvArray and build the recvList class
01671     for(int i(1), procID, size; i< recvArray[0]; i += size)
01672     {
01673       procID = recvArray[i++];
01674       size = recvArray[i++];
01675       recvList.push_back(exchange(procID, &recvArray[i], &recvArray[i+size], fieldDim) );
01676     }
01677 
01678     int me;
01679     MPI_Comm_rank ( comm, &me );
01680 
01681     exchangeList_Type::iterator it;
01682     for(it = recvList.begin(); it != recvList.end(); ++it )
01683     {
01684       if(it->procID == me) continue;
01685       MPI_Irecv(&(it->buffer[0]), it->buffer.size(), MPI_INT, it->procID, it->procID, comm, &it->reqID);
01686     }
01687 
01688     for(it = sendList.begin(); it != sendList.end(); ++it )
01689     {
01690       if(it->procID == me) continue;
01691       for(ID i=0; i< it->vec.size(); i++)
01692           for(int iComp=0; iComp<fieldDim; iComp++)
01693               it->buffer[fieldDim*i+iComp] = field[ fieldDim*it->vec[i] +iComp ];
01694 
01695       MPI_Isend(&(it->buffer[0]), it->buffer.size(), MPI_INT, it->procID, me, comm, &it->reqID);
01696     }
01697 
01698     for(it = recvList.begin(); it != recvList.end(); ++it )
01699     {
01700       if(it->procID == me) continue;
01701       MPI_Wait(&it->reqID, MPI_STATUS_IGNORE);
01702 
01703       for(int i=0; i< int(it->vec.size()); i++)
01704           for(int iComp=0; iComp<fieldDim; iComp++)
01705               field[ fieldDim*it->vec[i] + iComp] = it->buffer[fieldDim*i+iComp];
01706     }
01707 
01708     for(it = sendList.begin(); it != sendList.end(); ++it )
01709     {
01710       if(it->procID == me) continue;
01711       MPI_Wait(&it->reqID, MPI_STATUS_IGNORE);
01712     }
01713   }
01714 
01715   void allToAll(std::vector<int>& field, exchangeList_Type const * sendList, exchangeList_Type const * recvList, int fieldDim)
01716     {
01717         int me;
01718         MPI_Comm_rank ( comm, &me );
01719 
01720 
01721         for (int iComp = 0; iComp< fieldDim; iComp++)
01722         {
01723             exchangeList_Type::const_iterator it;
01724             for(it = recvList->begin(); it != recvList->end(); ++it )
01725             {
01726                 if(it->procID == me) continue;
01727                 MPI_Irecv(&(it->buffer[0]), it->buffer.size(), MPI_INT, it->procID, it->procID, comm, &it->reqID);
01728             }
01729 
01730             for(it = sendList->begin(); it != sendList->end(); ++it )
01731             {
01732                 if(it->procID == me) continue;
01733                 for(ID i=0; i< it->vec.size(); i++)
01734                     it->buffer[i] = field[ fieldDim*it->vec[i] +iComp ];
01735 
01736                 MPI_Isend(&(it->buffer[0]), it->buffer.size(), MPI_INT, it->procID, me, comm, &it->reqID);
01737             }
01738 
01739             for(it = recvList->begin(); it != recvList->end(); ++it )
01740             {
01741                 if(it->procID == me) continue;
01742                 MPI_Wait(&it->reqID, MPI_STATUS_IGNORE);
01743 
01744                 for(int i=0; i< int(it->vec.size()); i++)
01745                         field[ fieldDim*it->vec[i] + iComp] = it->buffer[i];
01746             }
01747 
01748             for(it = sendList->begin(); it != sendList->end(); ++it )
01749             {
01750                 if(it->procID == me) continue;
01751                 MPI_Wait(&it->reqID, MPI_STATUS_IGNORE);
01752             }
01753         }
01754     }
01755 
01756   void allToAll(double* field, exchangeList_Type const * sendList, exchangeList_Type const * recvList, int fieldDim)
01757     {
01758         int me;
01759         MPI_Comm_rank ( comm, &me );
01760 
01761         for (int iComp = 0; iComp< fieldDim; iComp++)
01762         {
01763             exchangeList_Type::const_iterator it;
01764   
01765             for(it = recvList->begin(); it != recvList->end(); ++it )
01766             {
01767                 if(it->procID == me) continue;
01768                 MPI_Irecv(&(it->doubleBuffer[0]), it->doubleBuffer.size(), MPI_DOUBLE, it->procID, it->procID, comm, &it->reqID);
01769             }
01770 
01771             for(it = sendList->begin(); it != sendList->end(); ++it )
01772             {
01773                 if(it->procID == me) continue;
01774                 for(ID i=0; i< it->vec.size(); i++)
01775                     it->doubleBuffer[i] = field[ fieldDim*it->vec[i] +iComp ];
01776 
01777                 MPI_Isend(&(it->doubleBuffer[0]), it->doubleBuffer.size(), MPI_DOUBLE, it->procID, me, comm, &it->reqID);
01778             }
01779 
01780             for(it = recvList->begin(); it != recvList->end(); ++it )
01781             {
01782                 if(it->procID == me) continue;
01783                 MPI_Wait(&it->reqID, MPI_STATUS_IGNORE);
01784 
01785                 for(int i=0; i< int(it->vec.size()); i++)
01786                         field[ fieldDim*it->vec[i] + iComp] = it->doubleBuffer[i];
01787             }
01788 
01789             for(it = sendList->begin(); it != sendList->end(); ++it )
01790             {
01791                 if(it->procID == me) continue;
01792                 MPI_Wait(&it->reqID, MPI_STATUS_IGNORE);
01793             }
01794         }
01795     }
01796 
01797   int initialize_iceProblem(int nTriangles)
01798     {
01799       bool keep_proc = nTriangles > 0;
01800 
01801       createReducedMPI(keep_proc, reducedComm);
01802 
01803     //    delete iceProblemPtr;
01804 
01805         isDomainEmpty = !keep_proc;
01806 
01807         // initialize ice problem pointer
01808         if(keep_proc)
01809         {
01810      //       iceProblemPtr = new ICEProblem(reducedComm);
01811             std::cout << nTriangles << " elements of the triangular grid are stored on this processor" <<std::endl;
01812         }
01813         else
01814         {
01815      //      iceProblemPtr = 0;
01816            std::cout << "No elements of the triangular grid are stored on this processor" <<std::endl;
01817         }
01818 
01819         return 0;
01820     }