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