| Directory: | ./ |
|---|---|
| File: | DegreeOfFreedom/dofBoundary.cpp |
| Date: | 2024-04-14 07:32:34 |
| Exec | Total | Coverage | |
|---|---|---|---|
| Lines: | 193 | 487 | 39.6% |
| Branches: | 206 | 1400 | 14.7% |
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // ______ ______ _ _ _____ ______ | ||
| 2 | // | ____| ____| | (_)/ ____| | ____| | ||
| 3 | // | |__ | |__ | | _| (___ ___| |__ | ||
| 4 | // | __| | __| | | | |\___ \ / __| __| | ||
| 5 | // | | | |____| |____| |____) | (__| |____ | ||
| 6 | // |_| |______|______|_|_____/ \___|______| | ||
| 7 | // Finite Elements for Life Sciences and Engineering | ||
| 8 | // | ||
| 9 | // License: LGL2.1 License | ||
| 10 | // FELiScE default license: LICENSE in root folder | ||
| 11 | // | ||
| 12 | // Main authors: J. Foulon & J-F. Gerbeau | ||
| 13 | // | ||
| 14 | |||
| 15 | // System includes | ||
| 16 | #include <unordered_map> | ||
| 17 | |||
| 18 | // External includes | ||
| 19 | |||
| 20 | // Project includes | ||
| 21 | #include "DegreeOfFreedom/dofBoundary.hpp" | ||
| 22 | #include "Solver/linearProblem.hpp" | ||
| 23 | #include "PETScInterface/petscMatrix.hpp" | ||
| 24 | |||
| 25 | namespace felisce | ||
| 26 | { | ||
| 27 | /*! \brief Class constructor. | ||
| 28 | It sets all flags to false | ||
| 29 | */ | ||
| 30 | 538 | DofBoundary::DofBoundary(): | |
| 31 | 538 | m_boundaryPatternAlreadyBuilt(false), | |
| 32 | 538 | m_dofBoundaryHasBeenInitialized(false), | |
| 33 | 538 | m_volumeInterfaceMappingBuilt(false) | |
| 34 | 538 | {} | |
| 35 | |||
| 36 | /*! \brief Initializer of the entire class | ||
| 37 | |||
| 38 | It has to be called in linear problem. | ||
| 39 | |||
| 40 | \param[in] lpb_ptr a pointer to the linear problem class. | ||
| 41 | |||
| 42 | The flag #m_dofBoundaryHasBeenInitialized is std::set to false. | ||
| 43 | */ | ||
| 44 | 8 | void DofBoundary::initialize(LinearProblem* lpb_ptr) { | |
| 45 | 8 | m_dofPtr = &lpb_ptr->dof(); | |
| 46 | 8 | m_ao = lpb_ptr->ao(); | |
| 47 | 8 | m_dofPartPtr = &lpb_ptr->dofPartition(); | |
| 48 | 8 | m_dofBoundaryHasBeenInitialized=true; | |
| 49 | 8 | } | |
| 50 | |||
| 51 | /*! \brief It creates the pattern for a sparse matrix defined on the boundary. | ||
| 52 | |||
| 53 | The pattern is taken directly from the volume pattern. | ||
| 54 | It implies: | ||
| 55 | - #m_dofBoundaryHasBeenInitialized = true | ||
| 56 | - #m_volumeInterfaceMappingBuilt = true. | ||
| 57 | In case it had already been called before it does not do anyting. | ||
| 58 | It raises the flag #m_boundaryPatternAlreadyBuilt to true. | ||
| 59 | |||
| 60 | #m_boundaryPattern is created. | ||
| 61 | */ | ||
| 62 | ✗ | void DofBoundary::buildBoundaryPattern() { | |
| 63 | ✗ | if ( m_boundaryPatternAlreadyBuilt ) | |
| 64 | ✗ | return; | |
| 65 | ✗ | if ( !m_dofBoundaryHasBeenInitialized ) | |
| 66 | ✗ | FEL_ERROR("dof boundary not yet initialized"); | |
| 67 | ✗ | if ( !m_volumeInterfaceMappingBuilt) | |
| 68 | ✗ | FEL_ERROR("volume interface mapping not yet build"); | |
| 69 | // if the pattern is already built than exit | ||
| 70 | ✗ | m_boundaryPatternAlreadyBuilt = true; | |
| 71 | |||
| 72 | // create a mapping to get local numbering of a dof given in global application numbering | ||
| 73 | // TODO check if the mapping is present in linearProblem | ||
| 74 | ✗ | felInt idLocal = 0; | |
| 75 | ✗ | std::unordered_map<felInt,felInt> glob2loc; | |
| 76 | ✗ | for (felInt i = 0; i < m_dofPtr->numDof(); i++) { | |
| 77 | ✗ | if ( (*m_dofPartPtr)[i] == MpiInfo::rankProc()) { | |
| 78 | ✗ | glob2loc[i] = idLocal; | |
| 79 | ✗ | idLocal++; | |
| 80 | } | ||
| 81 | } | ||
| 82 | |||
| 83 | // this structure will contain, for each node on the Boundary, the std::vector of the neighbours which belong to the boundary. | ||
| 84 | // it is local to the proc. | ||
| 85 | ✗ | std::vector< std::vector< felInt > > nodesNeighborhood( m_numLocalDofInterface ); | |
| 86 | // the number of nnz in the local pattern | ||
| 87 | ✗ | std::size_t nnz = 0; | |
| 88 | // building of nodesNeighborhood: | ||
| 89 | // for each surface node in this proc.. | ||
| 90 | ✗ | for ( int iNode = 0; iNode < m_numLocalDofInterface; iNode++ ) { | |
| 91 | // get all the neighbours! | ||
| 92 | // they can be both on surface and in the volume | ||
| 93 | // and in this proc or in another | ||
| 94 | ✗ | std::vector<felInt> currentVecOfNeighbour; | |
| 95 | ✗ | felInt iGlobAppli = m_loc2PetscVolume[iNode]; | |
| 96 | ✗ | AOPetscToApplication(m_ao,1,&iGlobAppli); | |
| 97 | ✗ | m_dofPtr->getNeighbourhood( glob2loc[iGlobAppli], currentVecOfNeighbour ); //currentVec is in global application numbering | |
| 98 | // we save only the neighbours that are on the surface | ||
| 99 | // we take also the neighbours on different proc | ||
| 100 | ✗ | for ( std::size_t iNeigh(0); iNeigh < currentVecOfNeighbour.size(); ++iNeigh ) { | |
| 101 | ✗ | felInt iPetsc = currentVecOfNeighbour[iNeigh]; | |
| 102 | ✗ | AOApplicationToPetsc(m_ao,1,&iPetsc); | |
| 103 | ✗ | if ( std::find( m_glob2PetscVolume.begin(), m_glob2PetscVolume.end(), iPetsc) != m_glob2PetscVolume.end() ) { | |
| 104 | ✗ | nodesNeighborhood[ iNode ].push_back(currentVecOfNeighbour[iNeigh]); | |
| 105 | } | ||
| 106 | } | ||
| 107 | // we update nnz | ||
| 108 | ✗ | nnz = nnz + nodesNeighborhood[ iNode ].size(); | |
| 109 | } | ||
| 110 | // resize of the pattern! | ||
| 111 | // nb of rows: m_numLocalDofInterface | ||
| 112 | // total number of non zeros | ||
| 113 | ✗ | m_boundaryPattern.resize(m_numLocalDofInterface,nnz); | |
| 114 | // we fill it! | ||
| 115 | ✗ | for ( int iNode = 0; iNode < m_numLocalDofInterface; iNode++) { | |
| 116 | ✗ | m_boundaryPattern.rowPointer(iNode+1) = m_boundaryPattern.rowPointer(iNode) + nodesNeighborhood[iNode].size(); | |
| 117 | ✗ | felInt pos = 0; | |
| 118 | ✗ | for(auto iCon = nodesNeighborhood[iNode].begin(); iCon != nodesNeighborhood[iNode].end(); iCon++) { | |
| 119 | // here we apply the std::unordered_map petscVolume2Glob to move the index *iCon from the global application volume ordering to | ||
| 120 | // the global application interface ordering | ||
| 121 | ✗ | felInt iPetscVol = *iCon; | |
| 122 | ✗ | AOApplicationToPetsc( m_ao, 1, &iPetscVol ); | |
| 123 | ✗ | m_boundaryPattern.columnIndex( m_boundaryPattern.rowPointer(iNode) + pos ) = m_petscVolume2Glob[ iPetscVol ]; | |
| 124 | ✗ | pos++; | |
| 125 | } | ||
| 126 | } | ||
| 127 | } | ||
| 128 | |||
| 129 | /*! \brief It allocates a sparse petsc matrix using the boundary pattern. | ||
| 130 | |||
| 131 | \param[in,out] theMatrix an empty matrix that will be allocated. | ||
| 132 | |||
| 133 | It is usefull, for instance, to allocate a mass matrix on the boarder. | ||
| 134 | It is necessary to pass a completely empty matrix. | ||
| 135 | |||
| 136 | If the boundary pattern was not built it builds it. | ||
| 137 | */ | ||
| 138 | ✗ | void DofBoundary::allocateMatrixOnBoundary( PetscMatrix& theMatrix ) { | |
| 139 | ✗ | if ( !m_boundaryPatternAlreadyBuilt ) | |
| 140 | ✗ | this->buildBoundaryPattern(); | |
| 141 | // only processors who actually own some nodes on the boundary | ||
| 142 | ✗ | if ( m_numLocalDofInterface > 0 ) { | |
| 143 | // we get the nb of procs that share the boundary nodes | ||
| 144 | felInt size; | ||
| 145 | ✗ | MPI_Comm_size( m_boundaryComm, &size ); | |
| 146 | |||
| 147 | // the matrix is now allocate in the same way as it is done for the system matrix | ||
| 148 | // parallel case | ||
| 149 | ✗ | if ( size > 1 ) { | |
| 150 | |||
| 151 | ✗ | std::vector<felInt> jCSR = m_boundaryPattern.columnIndices(); | |
| 152 | ✗ | AOApplicationToPetsc(m_aoInterface, m_boundaryPattern.numNonzeros(), jCSR.data()); | |
| 153 | |||
| 154 | ✗ | theMatrix.createAIJ(m_boundaryComm, m_numLocalDofInterface, m_numLocalDofInterface, | |
| 155 | m_numGlobalDofInterface, m_numGlobalDofInterface, 0, | ||
| 156 | FELISCE_PETSC_NULLPTR, 0, FELISCE_PETSC_NULLPTR); | ||
| 157 | ✗ | theMatrix.mpiAIJSetPreallocationCSR(m_boundaryPattern.rowPointer().data(), jCSR.data(), nullptr); | |
| 158 | ✗ | theMatrix.setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); | |
| 159 | ✗ | theMatrix.setFromOptions(); | |
| 160 | ✗ | } // serial case | |
| 161 | else { | ||
| 162 | |||
| 163 | ✗ | std::size_t numRows = m_boundaryPattern.numRows(); | |
| 164 | ✗ | std::vector<felInt> nnz( numRows ); | |
| 165 | ✗ | for ( std::size_t idof = 0; idof < numRows; idof++ ) { | |
| 166 | ✗ | nnz[idof] = m_boundaryPattern.numNonzerosInRow(idof); | |
| 167 | } | ||
| 168 | |||
| 169 | ✗ | theMatrix.createSeqAIJ(m_boundaryComm, numRows, numRows, 0, nnz.data()); | |
| 170 | ✗ | theMatrix.setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); | |
| 171 | ✗ | theMatrix.setFromOptions(); | |
| 172 | } | ||
| 173 | } | ||
| 174 | } | ||
| 175 | |||
| 176 | /***********************************************************************************/ | ||
| 177 | /***********************************************************************************/ | ||
| 178 | |||
| 179 | 308 | PetscVector DofBoundary::allocateBoundaryVector(vectorType type) | |
| 180 | { | ||
| 181 | // To use this function you have to firt initialize correctly dofBoundary | ||
| 182 | // look for instance into linearProblemNSRS.cpp | ||
| 183 | 308 | PetscVector v; | |
| 184 |
1/2✓ Branch 1 taken 308 times.
✗ Branch 2 not taken.
|
308 | if ( hasDofsOnBoundary() ) { |
| 185 | int numProcOnBoundary; | ||
| 186 |
1/2✓ Branch 2 taken 308 times.
✗ Branch 3 not taken.
|
308 | MPI_Comm_size( comm(), &numProcOnBoundary); |
| 187 |
2/4✓ Branch 0 taken 308 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 308 times.
|
308 | if ( numProcOnBoundary == 1 || type == sequential ) { |
| 188 | ✗ | v.createSeq( comm(), numGlobalDofInterface() ); | |
| 189 | } else { | ||
| 190 |
1/2✓ Branch 4 taken 308 times.
✗ Branch 5 not taken.
|
308 | v.createMPI( comm(), numLocalDofInterface(), numGlobalDofInterface() ); |
| 191 | } | ||
| 192 |
1/2✓ Branch 1 taken 308 times.
✗ Branch 2 not taken.
|
308 | v.setFromOptions(); |
| 193 | } | ||
| 194 | 308 | return v; | |
| 195 | } | ||
| 196 | |||
| 197 | /***********************************************************************************/ | ||
| 198 | /***********************************************************************************/ | ||
| 199 | |||
| 200 | 160 | void DofBoundary::restrictOnBoundary(PetscVector& volumeVector, PetscVector& parallelBoundaryVector) | |
| 201 | { | ||
| 202 |
1/2✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
|
160 | if ( hasDofsOnBoundary() ) { |
| 203 |
1/2✓ Branch 3 taken 160 times.
✗ Branch 4 not taken.
|
160 | std::vector<double> tmpValuesBD( numLocalDofInterface(), 0.0); |
| 204 |
1/2✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
|
160 | volumeVector.getValues(numLocalDofInterface(), loc2PetscVolPtr(), tmpValuesBD.data()); |
| 205 |
1/2✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
|
160 | parallelBoundaryVector.setValues(numLocalDofInterface(), loc2PetscBDPtr(), tmpValuesBD.data(), INSERT_VALUES); |
| 206 |
1/2✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
|
160 | parallelBoundaryVector.assembly(); |
| 207 | 160 | } | |
| 208 | 160 | } | |
| 209 | |||
| 210 | /***********************************************************************************/ | ||
| 211 | /***********************************************************************************/ | ||
| 212 | |||
| 213 | 152 | void DofBoundary::extendOnVolume(PetscVector& parallelVolumeVector, PetscVector& boundaryVector) | |
| 214 | { | ||
| 215 |
1/2✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
|
152 | if ( hasDofsOnBoundary() ) { |
| 216 |
1/2✓ Branch 3 taken 152 times.
✗ Branch 4 not taken.
|
152 | std::vector<double> tmpValuesBD( numLocalDofInterface(), 0.0); |
| 217 |
1/2✓ Branch 4 taken 152 times.
✗ Branch 5 not taken.
|
152 | boundaryVector.getValues( numLocalDofInterface(), loc2PetscBDPtr(), tmpValuesBD.data()); // Values are read from the boundary std::vector. |
| 218 |
1/2✓ Branch 4 taken 152 times.
✗ Branch 5 not taken.
|
152 | parallelVolumeVector.setValues( numLocalDofInterface(), loc2PetscVolPtr(), tmpValuesBD.data(), INSERT_VALUES); // Values are written in a PARALLEL volume std::vector. |
| 219 | 152 | } | |
| 220 | 152 | parallelVolumeVector.assembly(); | |
| 221 | 152 | } | |
| 222 | |||
| 223 | /***********************************************************************************/ | ||
| 224 | /***********************************************************************************/ | ||
| 225 | |||
| 226 | 16 | void DofBoundary::buildListOfBoundaryPetscDofs(LinearProblem* lpb_ptr, std::vector<int> labelOfInterface, int iUnknown, int iComponent) | |
| 227 | { | ||
| 228 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
16 | if ( ! m_dofBoundaryHasBeenInitialized ) |
| 229 | ✗ | FEL_ERROR("dof boundary not yet initialized"); | |
| 230 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | std::size_t iUC(m_dofPtr->getNumGlobComp(iUnknown,iComponent)); |
| 231 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | if ( iUC >= m_dofsHaveBeenComputed.size() ) { |
| 232 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | m_dofsHaveBeenComputed.resize( iUC+1, false); |
| 233 | ✗ | } else if ( m_dofsHaveBeenComputed[iUC] ) { | |
| 234 | ✗ | std::stringstream msg; | |
| 235 | ✗ | msg<<"list of boundary petsc dofs already computed for iUnknown: "<<iUnknown<<" iComp: "<<iComponent<<std::endl; | |
| 236 | ✗ | FEL_ERROR(msg.str().c_str()); | |
| 237 | } | ||
| 238 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | m_dofsHaveBeenComputed[iUC] = true; |
| 239 | |||
| 240 | 16 | const int idVar = lpb_ptr->listUnknown().idVariable(iUnknown); | |
| 241 | 16 | const int iMesh = lpb_ptr->listVariable()[idVar].idMesh(); | |
| 242 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16 times.
|
16 | if (FelisceParam::verbose() > 2 ) |
| 243 | ✗ | std::cout<<"["<<MpiInfo::rankProc()<<"]starting std::unordered_map interface creation in "<<lpb_ptr->name()<<" unknown id: "<<iUnknown<<", component: "<<iComponent<<std::endl; | |
| 244 | |||
| 245 | 16 | std::vector<setOfPairOfIntAndPoint> globalVecOfSet; | |
| 246 | |||
| 247 | // some useful std::vector | ||
| 248 | 16 | std::vector<Point*> elemPoint; | |
| 249 | 16 | std::vector<felInt> elemIdPoint; | |
| 250 | 16 | std::vector<felInt> vectorIdSupport; | |
| 251 | |||
| 252 | //to improve readability | ||
| 253 | 16 | GeometricMeshRegion::Pointer mesh = lpb_ptr->meshLocal(iMesh); | |
| 254 | 16 | const std::vector<GeometricMeshRegion::ElementType>& bagElementTypeDomainBoundary = mesh->bagElementTypeDomainBoundary(); | |
| 255 | |||
| 256 | //initialization of the counter numElement | ||
| 257 | felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ]; //std::vector of size 23 number of the different type of elements why we do not use a std::unordered_map | ||
| 258 |
2/2✓ Branch 0 taken 400 times.
✓ Branch 1 taken 16 times.
|
416 | for (int ityp=0; ityp < GeometricMeshRegion::m_numTypesOfElement; ityp++ ) { |
| 259 | 400 | numElement[ityp]=0; | |
| 260 | } | ||
| 261 | |||
| 262 | //========loop on the element type of boundary (e.g. triangles..not a real loop) | ||
| 263 |
2/2✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
|
32 | for (std::size_t i = 0; i < bagElementTypeDomainBoundary.size(); ++i) { |
| 264 | 16 | GeometricMeshRegion::ElementType eltType = bagElementTypeDomainBoundary[i]; | |
| 265 | |||
| 266 | 16 | int typeOfFiniteElement = lpb_ptr->listVariable()[idVar].finiteElementType(); | |
| 267 | 16 | const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second; | |
| 268 |
1/2✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
16 | const RefElement *refEle = geoEle->defineFiniteEle(eltType, typeOfFiniteElement, *mesh); |
| 269 |
2/4✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
|
16 | CurvilinearFiniteElement* fe = new CurvilinearFiniteElement(*refEle, *geoEle, lpb_ptr->listVariable()[idVar].degreeOfExactness()); |
| 270 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
|
16 | FEL_ASSERT(fe); |
| 271 |
3/6✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 16 times.
|
16 | if( FelisceParam::instance(lpb_ptr->instanceIndex()).flipNormal ) |
| 272 | ✗ | fe->m_sign = -1; | |
| 273 | |||
| 274 | 16 | const int numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType]; | |
| 275 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | elemPoint.resize(numPointPerElt, nullptr); |
| 276 |
1/2✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
|
16 | elemIdPoint.resize(numPointPerElt, 0); |
| 277 | |||
| 278 | //=======loop on subregion of the surface, with different labels. | ||
| 279 |
2/2✓ Branch 6 taken 40 times.
✓ Branch 7 taken 16 times.
|
56 | for(auto itRef = mesh->intRefToBegEndMaps[eltType].begin(); itRef != mesh->intRefToBegEndMaps[eltType].end(); itRef++) { |
| 280 | 40 | const int currentLabel = itRef->first; | |
| 281 | 40 | const int numEltPerLabel = itRef->second.second; | |
| 282 | |||
| 283 | 40 | bool Found( false ); //flag to select the correct label | |
| 284 |
5/6✓ Branch 4 taken 40 times.
✓ Branch 5 taken 40 times.
✓ Branch 6 taken 40 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 40 times.
✓ Branch 9 taken 40 times.
|
80 | for(auto it_labelNumber = labelOfInterface.begin(); it_labelNumber != labelOfInterface.end() && !Found; ++it_labelNumber) { |
| 285 |
2/2✓ Branch 1 taken 16 times.
✓ Branch 2 taken 24 times.
|
40 | if(*it_labelNumber == currentLabel) { |
| 286 | 16 | Found=true; | |
| 287 | |||
| 288 |
1/2✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
16 | setOfPairOfIntAndPoint localSetOfPointsAndIds(&Tools::lgPairComparison); |
| 289 | |||
| 290 |
2/2✓ Branch 0 taken 320 times.
✓ Branch 1 taken 16 times.
|
336 | for ( felInt iel = 0; iel < numEltPerLabel; iel++) { |
| 291 | // return each id of point of the element and coordinate in two arrays: | ||
| 292 | // elemPoint and elemIdPoint. | ||
| 293 |
1/2✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
|
320 | lpb_ptr->setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, vectorIdSupport); |
| 294 | |||
| 295 | // loop over all the support elements | ||
| 296 |
2/2✓ Branch 1 taken 320 times.
✓ Branch 2 taken 320 times.
|
640 | for (std::size_t it = 0; it < vectorIdSupport.size(); it++) { |
| 297 | // get the id of the support | ||
| 298 | 320 | felInt ielSupportDof = vectorIdSupport[it]; | |
| 299 |
2/2✓ Branch 1 taken 640 times.
✓ Branch 2 taken 320 times.
|
960 | for ( std::size_t iSurfDof(0); iSurfDof < (std::size_t) fe->numDof(); iSurfDof++ ) { |
| 300 | |||
| 301 | felInt iDofGlobScalar; | ||
| 302 |
1/2✓ Branch 2 taken 640 times.
✗ Branch 3 not taken.
|
640 | lpb_ptr->dof().loc2glob( ielSupportDof, iSurfDof, idVar, iComponent, iDofGlobScalar); |
| 303 |
1/2✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
|
640 | AOApplicationToPetsc(m_ao,1,&iDofGlobScalar); |
| 304 | Point P; | ||
| 305 |
1/2✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
|
640 | lpb_ptr->getSupportCoordinate(eltType, numElement[eltType], iSurfDof, iUnknown, P); |
| 306 |
2/4✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 640 times.
✗ Branch 5 not taken.
|
640 | localSetOfPointsAndIds.insert(std::make_pair(iDofGlobScalar,P) ); |
| 307 | } | ||
| 308 | } | ||
| 309 | 320 | numElement[eltType]++; | |
| 310 | } | ||
| 311 | |||
| 312 |
2/2✓ Branch 1 taken 8 times.
✓ Branch 2 taken 8 times.
|
16 | if ( int(m_localVecOfSet.size()) <= currentLabel ) { |
| 313 |
2/4✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
|
8 | m_localVecOfSet.resize(currentLabel+1, setOfPairOfIntAndPoint(&Tools::lgPairComparison) ); |
| 314 | } | ||
| 315 |
1/2✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
16 | m_localVecOfSet[currentLabel] = localSetOfPointsAndIds; |
| 316 | 16 | } | |
| 317 | } | ||
| 318 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 16 times.
|
40 | if (!Found) |
| 319 | 24 | numElement[eltType] += numEltPerLabel; | |
| 320 | } | ||
| 321 | |||
| 322 |
1/2✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
|
16 | delete fe; |
| 323 | }//end of the loop | ||
| 324 | |||
| 325 |
1/2✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
16 | const int idUnknComp = lpb_ptr->dof().getNumGlobComp(iUnknown,iComponent); |
| 326 |
2/2✓ Branch 4 taken 16 times.
✓ Branch 5 taken 16 times.
|
32 | for(auto it_labelNumber = labelOfInterface.begin(); it_labelNumber != labelOfInterface.end(); ++it_labelNumber) { |
| 327 | |||
| 328 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16 times.
|
16 | if ( FelisceParam::verbose() > 3 ) |
| 329 | ✗ | std::cout<<"["<<MpiInfo::rankProc()<<"] starting gathering phase"<<std::endl; | |
| 330 | |||
| 331 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
|
16 | if ( int( m_localVecOfSet.size() ) <= *it_labelNumber ) |
| 332 | ✗ | m_localVecOfSet.resize(*it_labelNumber+1, setOfPairOfIntAndPoint(&Tools::lgPairComparison) ); | |
| 333 | |||
| 334 |
1/2✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
|
16 | if ( int( globalVecOfSet.size() ) <= *it_labelNumber ) |
| 335 |
2/4✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
|
16 | globalVecOfSet.resize(*it_labelNumber+1, setOfPairOfIntAndPoint(&Tools::lgPairComparison) ); |
| 336 | |||
| 337 |
1/2✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
|
16 | Tools::allGatherSetPair(m_localVecOfSet[*it_labelNumber],globalVecOfSet[*it_labelNumber]); |
| 338 | 16 | const auto& r_global_vec = globalVecOfSet[*it_labelNumber]; | |
| 339 | |||
| 340 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16 times.
|
16 | if ( FelisceParam::verbose() > 3 ) |
| 341 | ✗ | std::cout<<"["<<MpiInfo::rankProc()<<"] ending gathering phase"<<std::endl; | |
| 342 | |||
| 343 | 16 | int cPoint=0; | |
| 344 | |||
| 345 |
2/2✓ Branch 3 taken 1296 times.
✓ Branch 4 taken 16 times.
|
1312 | for(auto it= r_global_vec.begin(); it!= r_global_vec.end(); ++it) { |
| 346 |
3/6✓ Branch 2 taken 1296 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1296 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1296 times.
✗ Branch 10 not taken.
|
1296 | m_listOfBoundaryPetscDofs[ idUnknComp ][*it_labelNumber][cPoint] = it->first; |
| 347 | 1296 | cPoint++; | |
| 348 | } | ||
| 349 | 16 | cPoint=0; | |
| 350 |
3/4✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 8 times.
|
16 | if ( m_mapLab2AllPoints[*it_labelNumber].size() == 0 ) { |
| 351 |
2/2✓ Branch 4 taken 648 times.
✓ Branch 5 taken 8 times.
|
656 | for(auto it= r_global_vec.begin(); it!= r_global_vec.end(); ++it) { |
| 352 |
1/2✓ Branch 0 taken 648 times.
✗ Branch 1 not taken.
|
648 | if (cPoint==0) { |
| 353 |
1/2✓ Branch 3 taken 648 times.
✗ Branch 4 not taken.
|
648 | m_mapLab2FirstPoint[*it_labelNumber] = it->second; |
| 354 | } | ||
| 355 |
2/4✓ Branch 2 taken 648 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 648 times.
✗ Branch 7 not taken.
|
648 | m_mapLab2AllPoints[*it_labelNumber].push_back(it->second); |
| 356 | } | ||
| 357 | 8 | cPoint++; | |
| 358 | } | ||
| 359 | |||
| 360 | } | ||
| 361 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16 times.
|
16 | if ( FelisceParam::verbose() > 3 ) |
| 362 | ✗ | std::cout<<"["<<MpiInfo::rankProc()<<"] ended"<<std::endl; | |
| 363 |
2/4✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16 times.
|
16 | if ( FelisceParam::verbose() > 3 ) |
| 364 | ✗ | std::cout<<"["<<MpiInfo::rankProc()<<"] "<<m_listOfBoundaryPetscDofs[idUnknComp].size()<<" labels, for "<<idUnknComp<<" id unknown and component"<<std::endl; | |
| 365 | 16 | } | |
| 366 | |||
| 367 | /***********************************************************************************/ | ||
| 368 | /***********************************************************************************/ | ||
| 369 | |||
| 370 | ✗ | void DofBoundary::writeSurfaceInterpolatorFile(int iUnknown, int numComp, std::string filename, std::string folder) const { | |
| 371 | // iUnknwon: index of the unknown we want to interpolate. | ||
| 372 | // numComp: number of components to distinguish between scalars (1) and vectors (2 or 3). | ||
| 373 | |||
| 374 | // Goal of the function is to write a list of the points, belonging to the surface, | ||
| 375 | // that support a degree of freedom (or more) for this unknown. | ||
| 376 | |||
| 377 | // ids: for each point we store the id (in globalApplicationOrderingBD) or the ids in the case of a std::vector(that's why a std::set is used.) | ||
| 378 | ✗ | std::vector<std::vector<int> > ids; | |
| 379 | // labels: for each point we save the surfaceLabel. | ||
| 380 | // if he belongs to more than one label we store the first one found. | ||
| 381 | ✗ | std::vector<felInt> labels; | |
| 382 | // points: for each point we actually store the Point. | ||
| 383 | ✗ | std::vector<Point> points; | |
| 384 | |||
| 385 | // Checking that everything was allocated | ||
| 386 | ✗ | for ( felInt iComp(0); iComp<numComp; ++iComp) { | |
| 387 | ✗ | std::size_t iUC=m_dofPtr->getNumGlobComp(iUnknown, iComp); | |
| 388 | ✗ | if ( iUC >= m_dofsHaveBeenComputed.size() || !m_dofsHaveBeenComputed[iUC] ) { | |
| 389 | ✗ | std::stringstream msg; | |
| 390 | ✗ | msg<<"Error in writing surface interpolator file for iUnknown: "<<iUnknown<<" and iComp: "<<iComp<<std::endl; | |
| 391 | ✗ | FEL_ERROR(msg.str().c_str()); | |
| 392 | } | ||
| 393 | } | ||
| 394 | |||
| 395 | // We do a first loop on all the points considering only the first component. | ||
| 396 | // The idea is that we need to build a sort of index for the points and it is not easy to do it | ||
| 397 | // since the point class does not work very well as key-value for a std::unordered_map. | ||
| 398 | |||
| 399 | ✗ | felInt iUC0 = m_dofPtr->getNumGlobComp(iUnknown, /*iComp*/0); | |
| 400 | |||
| 401 | // The idea is that for all the components we will loop over the labels and their points in the same way. | ||
| 402 | // When looping for the first component we associate an incremental index (cPointUnique) to the points | ||
| 403 | // that we find. If passing to a different label we find the same point twice we are able to tell it since, | ||
| 404 | // for the same component the point will be associated to the same degreeOfFreedom. | ||
| 405 | // The list of already found deegreesOfFreedom is saved in a std::set (alreadyFound). | ||
| 406 | // For each point we come across either we put the cPointUnique index in the std::vector idForiUC0 | ||
| 407 | // or we put (-1) if the point was already found in a different label. | ||
| 408 | ✗ | int cPointUnique(0); | |
| 409 | ✗ | std::set<int> alreadyFound; | |
| 410 | ✗ | std::vector<int> idForiUC0; | |
| 411 | |||
| 412 | // For each label | ||
| 413 | ✗ | for ( auto itLabel = m_listOfBoundaryPetscDofs.at(iUC0).begin(); itLabel != m_listOfBoundaryPetscDofs.at(iUC0).end(); ++itLabel) { | |
| 414 | // Point index label-based | ||
| 415 | ✗ | int cpoint(0); | |
| 416 | // for each point | ||
| 417 | ✗ | for ( auto itDof = itLabel->second.begin(); itDof != itLabel->second.end(); ++itDof) { | |
| 418 | ✗ | FEL_ASSERT(cpoint==itDof->first); | |
| 419 | // We insert the dof id in the list | ||
| 420 | ✗ | auto ret = alreadyFound.insert(m_petscVolume2Glob.at(itDof->second)); | |
| 421 | // We check if we already found this dof. | ||
| 422 | ✗ | if ( !ret.second ) { | |
| 423 | // In this case we do nothing | ||
| 424 | // and we add a minus one to idForiUC0 | ||
| 425 | ✗ | idForiUC0.push_back(-1); | |
| 426 | } else { | ||
| 427 | // In case of a new point | ||
| 428 | // we add the point index | ||
| 429 | ✗ | idForiUC0.push_back(cPointUnique); | |
| 430 | ✗ | cPointUnique++; | |
| 431 | // we save the point | ||
| 432 | ✗ | points.push_back(m_mapLab2AllPoints.find(itLabel->first)->second[cpoint]); | |
| 433 | // ...the label | ||
| 434 | ✗ | labels.push_back(itLabel->first); | |
| 435 | // ...and the dof id | ||
| 436 | ✗ | std::vector<int> tmp; | |
| 437 | ✗ | tmp.push_back(m_petscVolume2Glob.at(itDof->second)); | |
| 438 | ✗ | ids.push_back(tmp); | |
| 439 | } | ||
| 440 | //increase the index label-based | ||
| 441 | ✗ | cpoint++; | |
| 442 | } | ||
| 443 | } | ||
| 444 | |||
| 445 | // Now we are ready for a second round: | ||
| 446 | // if there are more components to interpolate | ||
| 447 | // we loop again in the same way as we just did for the first component | ||
| 448 | |||
| 449 | ✗ | for ( felInt iComp(1); iComp<numComp; ++iComp) { | |
| 450 | // This counter will be used to acces idForiUC0 and either extract -1 | ||
| 451 | // or extract the cPointUnique index. | ||
| 452 | ✗ | int count(0); | |
| 453 | ✗ | std::size_t iUC=m_dofPtr->getNumGlobComp(iUnknown, iComp); | |
| 454 | // for each label | ||
| 455 | ✗ | for ( auto itLabel = m_listOfBoundaryPetscDofs.at(iUC).begin(); itLabel != m_listOfBoundaryPetscDofs.at(iUC).end(); ++itLabel) { | |
| 456 | // counter label-based | ||
| 457 | ✗ | int cpoint=0; | |
| 458 | // loop over the points of this label | ||
| 459 | ✗ | for ( auto itDof = itLabel->second.begin(); itDof != itLabel->second.end(); ++itDof) { | |
| 460 | // verify if we have to skip this point | ||
| 461 | ✗ | if ( idForiUC0[count] >= 0 ) { | |
| 462 | // insert the dof id of the current component | ||
| 463 | ✗ | ids[idForiUC0[count]].push_back(m_petscVolume2Glob.at(itDof->second)); | |
| 464 | ✗ | FEL_ASSERT(Tools::equalTol(points[idForiUC0[count]].x(),m_mapLab2AllPoints.find(itLabel->first)->second[cpoint].x(),1e-16)); | |
| 465 | ✗ | FEL_ASSERT(Tools::equalTol(points[idForiUC0[count]].y(),m_mapLab2AllPoints.find(itLabel->first)->second[cpoint].y(),1e-16)); | |
| 466 | ✗ | FEL_ASSERT(Tools::equalTol(points[idForiUC0[count]].z(),m_mapLab2AllPoints.find(itLabel->first)->second[cpoint].z(),1e-16)); | |
| 467 | } | ||
| 468 | ✗ | count++; | |
| 469 | ✗ | cpoint++; | |
| 470 | } | ||
| 471 | } | ||
| 472 | } | ||
| 473 | // Opening the file | ||
| 474 | ✗ | std::stringstream filenameFull,cmd; | |
| 475 | ✗ | cmd<<"mkdir -p "<<folder; | |
| 476 | ✗ | int ierr = system(cmd.str().c_str()); | |
| 477 | ✗ | if (ierr>0) FEL_ERROR("error in execution of " + cmd.str() ); | |
| 478 | ✗ | filenameFull<<folder <<"/"<< filename; | |
| 479 | ✗ | std::ofstream surfaceInterpolatorFile; | |
| 480 | ✗ | surfaceInterpolatorFile.open(filenameFull.str().c_str()); | |
| 481 | ✗ | if ( ! surfaceInterpolatorFile.is_open() ) FEL_ERROR("not able to open the file"); | |
| 482 | |||
| 483 | // Writing the file | ||
| 484 | ✗ | surfaceInterpolatorFile<<points.size()<<std::endl; | |
| 485 | ✗ | for ( std::size_t curPoint(0); curPoint<points.size(); ++curPoint) { | |
| 486 | ✗ | surfaceInterpolatorFile<<numComp<<" "<<points[curPoint].x()<<" "<<points[curPoint].y()<<" "<<points[curPoint].z()<<" "; | |
| 487 | ✗ | surfaceInterpolatorFile<<labels[curPoint]<<" "; | |
| 488 | ✗ | FEL_ASSERT((std::size_t)numComp==ids[curPoint].size()); | |
| 489 | ✗ | for(auto it=ids[curPoint].begin(); it!= ids[curPoint].end(); it++) { | |
| 490 | ✗ | surfaceInterpolatorFile<<*it<<" "; | |
| 491 | } | ||
| 492 | ✗ | surfaceInterpolatorFile<<std::endl; | |
| 493 | } | ||
| 494 | ✗ | surfaceInterpolatorFile.close(); | |
| 495 | } | ||
| 496 | |||
| 497 | /***********************************************************************************/ | ||
| 498 | /***********************************************************************************/ | ||
| 499 | |||
| 500 | 1 | void DofBoundary::exportInterfacePoints(int iUnknown, int iComp) const { | |
| 501 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | std::size_t iUC=m_dofPtr->getNumGlobComp(iUnknown, iComp); |
| 502 |
4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
|
1 | if ( iUC >= m_dofsHaveBeenComputed.size() || !m_dofsHaveBeenComputed[iUC] ) { |
| 503 | ✗ | std::stringstream msg; | |
| 504 | ✗ | msg<<"trying to export points for iUnknown: "<<iUnknown<<" and iComp: "<<iComp<<std::endl; | |
| 505 | ✗ | FEL_ERROR(msg.str().c_str()); | |
| 506 | } | ||
| 507 | /*! a matlab file "interface.m" is written into the result folder | ||
| 508 | to import it is enough to run such file in a matlab shell | ||
| 509 | it will create a matrix called "interfacePoints" with | ||
| 510 | three columns and a row for each point | ||
| 511 | the order of the point is the order of the application | ||
| 512 | */ | ||
| 513 | 1 | std::unordered_map<int,Point> id2Point; | |
| 514 |
4/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 1 times.
|
2 | for(auto itLabel = m_listOfBoundaryPetscDofs.at(iUC).begin(); itLabel != m_listOfBoundaryPetscDofs.at(iUC).end(); ++itLabel) { |
| 515 | 1 | int cpoint(0); | |
| 516 |
2/2✓ Branch 5 taken 81 times.
✓ Branch 6 taken 1 times.
|
82 | for(auto itDof = itLabel->second.begin(); itDof != itLabel->second.end(); ++itDof) { |
| 517 |
2/4✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 8 taken 81 times.
✗ Branch 9 not taken.
|
81 | id2Point[itDof->second]=m_mapLab2AllPoints.find(itLabel->first)->second[cpoint]; |
| 518 |
2/4✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
|
81 | if ( FelisceParam::verbose() > 2 ) |
| 519 | ✗ | std::cout<<" id: "<<itDof->second<<" point "<<m_mapLab2AllPoints.find(itLabel->first)->second[cpoint]<<std::endl; | |
| 520 | 81 | cpoint++; | |
| 521 | } | ||
| 522 | } | ||
| 523 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | std::stringstream filename,cmd; |
| 524 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
1 | cmd<<"mkdir -p "<<FelisceParam::instance().resultDir; |
| 525 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
1 | int ierr = system(cmd.str().c_str()); |
| 526 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (ierr>0) { |
| 527 | ✗ | FEL_ERROR("error in execution of " + cmd.str() ); | |
| 528 | } | ||
| 529 | |||
| 530 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
1 | filename<< FelisceParam::instance().resultDir << "interface.m"; |
| 531 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | std::ofstream interfaceFile; |
| 532 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
1 | interfaceFile.open(filename.str().c_str()); |
| 533 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
|
1 | if ( ! interfaceFile.is_open() ) |
| 534 | ✗ | FEL_ERROR("not able to open the file"); | |
| 535 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | interfaceFile<<"interfacePoints = [ "; |
| 536 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 1 times.
|
82 | for ( int i(0); i<m_numGlobalDofInterface; ++i) { |
| 537 |
2/4✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
|
81 | if ( id2Point.count(m_glob2PetscVolume[i]) > 0 ) { |
| 538 |
2/4✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 81 times.
|
81 | if ( FelisceParam::verbose() > 2 ) |
| 539 | ✗ | std::cout<<" id2: "<<m_glob2PetscVolume[i]<<" point "<<id2Point[m_glob2PetscVolume[i]]<<std::endl; | |
| 540 |
1/2✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
|
81 | Point P=id2Point[m_glob2PetscVolume[i]]; |
| 541 |
6/12✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 81 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 81 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 81 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 81 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 81 times.
✗ Branch 20 not taken.
|
81 | interfaceFile<<P.x()<<" , "<<P.y()<<" , "<<P.z()<<" ; "; |
| 542 | } | ||
| 543 | } | ||
| 544 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | interfaceFile<<"];"<<std::endl; |
| 545 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | interfaceFile.close(); |
| 546 | 1 | } | |
| 547 | |||
| 548 | // for compatibility with previous definition of buildBoundaryVolumeMapping, but also | ||
| 549 | // to have a better code in the linear problems: | ||
| 550 | 4 | void DofBoundary::buildBoundaryVolumeMapping( int iUnknown, int iComponent ) { | |
| 551 |
2/4✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
|
4 | this->buildBoundaryVolumeMapping( iUnknown, std::vector<int>(1,iComponent) ); |
| 552 | 4 | } | |
| 553 | |||
| 554 | // One mapping at the time is allowed | ||
| 555 | // it means that you choose an unknown and a (std::set of) component and you create a mapping! | ||
| 556 | // for instance in Navier Stokes you can build mapping between the ?whole? volume problem ( vel x 3 + pre ) and | ||
| 557 | // one component of the velocity OR the pressure on the interface. | ||
| 558 | 4 | void DofBoundary::buildBoundaryVolumeMapping( int iUnknown, std::vector<int> components ) { | |
| 559 | |||
| 560 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if ( m_volumeInterfaceMappingBuilt ) { |
| 561 | #ifndef NDEBUG | ||
| 562 | ✗ | std::stringstream msg; | |
| 563 | ✗ | msg<<"trying to call builBoundaryVolumeMapping again with iUnknwon: "<<iUnknown<<" and component(s): "; | |
| 564 | ✗ | for(auto itComponent=components.begin(); itComponent!=components.end(); ++itComponent) { | |
| 565 | ✗ | msg<<*itComponent<<" "; | |
| 566 | } | ||
| 567 | ✗ | msg<<std::endl; | |
| 568 | ✗ | FEL_WARNING(msg.str().c_str()); | |
| 569 | #endif | ||
| 570 | |||
| 571 | ✗ | return; | |
| 572 | } | ||
| 573 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
|
4 | if ( ! m_dofBoundaryHasBeenInitialized ) |
| 574 | ✗ | FEL_ERROR("dof boundary not yet initialized"); | |
| 575 | 4 | m_volumeInterfaceMappingBuilt=true; | |
| 576 | |||
| 577 | 4 | std::set<int> setGlobVolPetscIds; | |
| 578 |
2/2✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
|
8 | for (std::size_t idComponent(0); idComponent<components.size(); ++idComponent) { |
| 579 | 4 | int iComponent = components[idComponent]; | |
| 580 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | int iUnknComp = m_dofPtr->getNumGlobComp(iUnknown, iComponent); |
| 581 | |||
| 582 | // we check some basic assumptions | ||
| 583 |
2/8✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
4 | FEL_CHECK( m_listOfBoundaryPetscDofs[iUnknComp].size() > 0, "not even one label stored in m_listOfBoundaryPetscDofs" ); |
| 584 |
2/8✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
4 | FEL_CHECK( m_listOfBoundaryPetscDofs[iUnknComp].begin()->second.size() > 0, "the first label stored in the std::unordered_map is empty" ); |
| 585 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
|
4 | if ( FelisceParam::verbose() > 2) { |
| 586 | ✗ | PetscPrintf( MpiInfo::petscComm(), "Starting creation of all the mappings to define petsc objects in the interface\n"); | |
| 587 | } | ||
| 588 | |||
| 589 | // =========================================================== | ||
| 590 | // the global volume petsc ids are stored in the m_listOfBoundaryPetscDofs | ||
| 591 | // container depending on their label, therefore some of those | ||
| 592 | // ids could be repeted if they are at the interface | ||
| 593 | // between different labels | ||
| 594 | // we "gather" all the ids in one std::set to destroy | ||
| 595 | // duplicates and to sort them | ||
| 596 |
4/6✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 4 times.
✓ Branch 11 taken 4 times.
|
8 | for(auto itLabel = m_listOfBoundaryPetscDofs[iUnknComp].begin(); itLabel != m_listOfBoundaryPetscDofs[iUnknComp].end(); ++itLabel) { |
| 597 |
2/2✓ Branch 6 taken 324 times.
✓ Branch 7 taken 4 times.
|
328 | for(auto itDof = itLabel->second.begin(); itDof != itLabel->second.end(); ++itDof) { |
| 598 |
1/2✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
|
324 | setGlobVolPetscIds.insert(itDof->second); |
| 599 | } | ||
| 600 | } | ||
| 601 | } | ||
| 602 | // =========================================================== | ||
| 603 | // m_glob2PetscVolume = | ||
| 604 | // =========================================================== | ||
| 605 | // TODO: this is simply a copy of the std::set into a std::vector, there is for sure a better way of doing it. | ||
| 606 | // we allocate the std::unordered_map | ||
| 607 | // applicationInterface2PetscVol (m_glob2PetscVolume), actually it is a std::vector of integers and | ||
| 608 | // it maps: | ||
| 609 | // ( 0 , m_numGlobalDofInterface - 1) --> (???) \in global volume petsc ids | ||
| 610 | // so that from the global interface application id you can get the global volume petsc id | ||
| 611 |
2/2✓ Branch 4 taken 324 times.
✓ Branch 5 taken 4 times.
|
328 | for(auto itDof = setGlobVolPetscIds.begin(); itDof != setGlobVolPetscIds.end(); ++itDof) { |
| 612 |
1/2✓ Branch 2 taken 324 times.
✗ Branch 3 not taken.
|
324 | m_glob2PetscVolume.push_back(*itDof); // global petsc volume id = m_glob2PetscVolume[global application interface id] |
| 613 | } | ||
| 614 | 4 | setGlobVolPetscIds.clear(); | |
| 615 | 4 | m_numGlobalDofInterface = m_glob2PetscVolume.size(); | |
| 616 | |||
| 617 | |||
| 618 | // =========================================================== | ||
| 619 | // m_petscVolume2Glob = | ||
| 620 | // =========================================================== | ||
| 621 | 4 | int cpoint=0; | |
| 622 |
2/2✓ Branch 3 taken 324 times.
✓ Branch 4 taken 4 times.
|
328 | for(auto petscVolId=m_glob2PetscVolume.begin(); petscVolId != m_glob2PetscVolume.end(); ++petscVolId, ++cpoint) { |
| 623 | 324 | int aus( *petscVolId ); | |
| 624 |
1/2✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
|
324 | m_petscVolume2Glob[ aus ] = cpoint; |
| 625 | } | ||
| 626 | |||
| 627 | // ============================================================== | ||
| 628 | // m_loc2GlobInterface & m_loc2PetscVolume = | ||
| 629 | // ============================================================== | ||
| 630 | // Here we create the std::vector, | ||
| 631 | // local to the proc, that contains the global application | ||
| 632 | // ids of the dof at the interface, i.e, | ||
| 633 | // m_loc2GlobInterface: | ||
| 634 | // ( 0 , m_numLocalDofInterface - 1 ) --> ( 0 , m_numGlobalDofInterface - 1 ) | ||
| 635 | // we create an auxiliary std::vector initialized with the global | ||
| 636 | // volume petsc ids and then we use m_ao | ||
| 637 | // to replace those values with global colume application ids | ||
| 638 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | std::vector<felInt> applicationIdsVolume(m_glob2PetscVolume); |
| 639 |
1/2✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
4 | AOPetscToApplication(m_ao,m_numGlobalDofInterface,&applicationIdsVolume[0]); |
| 640 | // this is done because we want to check if those dofs belong | ||
| 641 | // to this proc or not and the std::vector m_dofPart | ||
| 642 | // is defined on the global volume application ids not the petsc ids | ||
| 643 | // at the same time we define a variable startCounter | ||
| 644 | // which is equal to the number of dofs in the processors | ||
| 645 | // with a lower rank than the current | ||
| 646 | // it is zero for the proc 0 | ||
| 647 | // and it should be numGlob-numLoc-1 for the last proc | ||
| 648 | 4 | felInt startCounter(0); | |
| 649 |
2/2✓ Branch 0 taken 324 times.
✓ Branch 1 taken 4 times.
|
328 | for ( felInt iGlobDof(0); iGlobDof < m_numGlobalDofInterface; ++iGlobDof ) { |
| 650 | 324 | felInt dofProc( (*m_dofPartPtr)[ applicationIdsVolume[iGlobDof] ]); | |
| 651 |
3/4✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 81 times.
✓ Branch 4 taken 243 times.
|
324 | if ( dofProc == MpiInfo::rankProc() ) { |
| 652 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | m_loc2GlobInterface.push_back(iGlobDof); |
| 653 | } | ||
| 654 |
3/4✓ Branch 1 taken 324 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 105 times.
✓ Branch 4 taken 219 times.
|
324 | if ( dofProc < MpiInfo::rankProc() ) { |
| 655 | 105 | ++startCounter; | |
| 656 | } | ||
| 657 | } | ||
| 658 | 4 | applicationIdsVolume.clear(); | |
| 659 | 4 | m_numLocalDofInterface = m_loc2GlobInterface.size(); | |
| 660 |
2/2✓ Branch 4 taken 81 times.
✓ Branch 5 taken 4 times.
|
85 | for(auto iLocDof = m_loc2GlobInterface.begin(); iLocDof != m_loc2GlobInterface.end(); ++iLocDof ) { |
| 661 |
1/2✓ Branch 3 taken 81 times.
✗ Branch 4 not taken.
|
81 | m_loc2PetscVolume.push_back(m_glob2PetscVolume[*iLocDof]); |
| 662 | } | ||
| 663 | |||
| 664 |
4/8✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
|
4 | MPI_Comm_split(MpiInfo::petscComm(), ( m_numLocalDofInterface > 0 ) ? 1 : 0, MpiInfo::rankProc(), &m_boundaryComm); |
| 665 | int rankB,sizeB; | ||
| 666 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | MPI_Comm_rank(m_boundaryComm, &rankB); |
| 667 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | MPI_Comm_size(m_boundaryComm, &sizeB); |
| 668 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
|
4 | if (FelisceParam::verbose()>2) |
| 669 | ✗ | std::cout<<"["<<MpiInfo::rankProc()<<"] "<<"rank on the boarder "<< rankB << " size of this group " << sizeB << std::endl; | |
| 670 | |||
| 671 | // ============================================================== | ||
| 672 | // m_aoInterface = | ||
| 673 | // ============================================================== | ||
| 674 | // we create the std::vector petscIds | ||
| 675 | // which contains ( startCounter, startCounter + 1, ..., startCounter + numLocDof -1 ) | ||
| 676 | 4 | std::vector<felInt> petscIds; | |
| 677 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
4 | petscIds.reserve(m_numLocalDofInterface); |
| 678 | |||
| 679 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 4 times.
|
85 | for ( felInt iLocDof(0); iLocDof<m_numLocalDofInterface; ++iLocDof ) { |
| 680 |
1/2✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
|
81 | petscIds.push_back(startCounter + iLocDof); |
| 681 | } | ||
| 682 | // we build m_aoInterface | ||
| 683 | // this will allow you to move from the global petsc ordering defined in petscIds and | ||
| 684 | // the global application ordering. | ||
| 685 | // consider that in serial this is the identity mapping | ||
| 686 | // but in parallel nodes are contiguous in petsc ordering | ||
| 687 |
1/2✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
|
4 | AOCreateBasic(m_boundaryComm, m_numLocalDofInterface, &m_loc2GlobInterface[0], &petscIds[0], &m_aoInterface); |
| 688 | 4 | petscIds.clear(); | |
| 689 | |||
| 690 | |||
| 691 | // ============================================================== | ||
| 692 | // m_loc2PetscInterface = | ||
| 693 | // ============================================================== | ||
| 694 | // from local index, directly to the petsc index of the interface | ||
| 695 |
2/2✓ Branch 0 taken 81 times.
✓ Branch 1 taken 4 times.
|
85 | for ( felInt iDofInterface(0); iDofInterface < m_numLocalDofInterface; ++iDofInterface) { |
| 696 |
1/2✓ Branch 2 taken 81 times.
✗ Branch 3 not taken.
|
81 | m_loc2PetscInterface.push_back(m_loc2GlobInterface[iDofInterface]); |
| 697 | } | ||
| 698 |
1/2✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
|
4 | AOApplicationToPetsc(m_aoInterface,m_numLocalDofInterface,&m_loc2PetscInterface[0]); |
| 699 | |||
| 700 | // bool m_exportInterface=true; | ||
| 701 | // if ( m_exportInterface ) { | ||
| 702 |
3/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 3 times.
|
4 | if ( MpiInfo::rankProc() == 0 ) { |
| 703 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | this->exportInterfacePoints( iUnknown, components[0] );//we export the points only for the first component. |
| 704 | } | ||
| 705 | // } | ||
| 706 |
2/4✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
|
4 | if ( FelisceParam::verbose() > 2 ) { |
| 707 | ✗ | displayBoundaryVolumeMapping(); | |
| 708 | } | ||
| 709 | 4 | } | |
| 710 | |||
| 711 | /*! \brief post-process function used in the examples | ||
| 712 | It writes a readable file to be plot in python. | ||
| 713 | */ | ||
| 714 | ✗ | void DofBoundary::BoundaryVolumeMappingForExample() const { | |
| 715 | ✗ | if ( ! m_volumeInterfaceMappingBuilt ) | |
| 716 | ✗ | FEL_ERROR("volume interface mapping not yet build"); | |
| 717 | |||
| 718 | //! The master processor start to write: | ||
| 719 | ✗ | if ( MpiInfo::rankProc() == 0 ) { | |
| 720 | //! a file BVMapping.dat is opened | ||
| 721 | ✗ | std::stringstream cmd; | |
| 722 | ✗ | cmd<<"mkdir -p "<<FelisceParam::instance().resultDir; | |
| 723 | ✗ | int ierr = system( cmd.str().c_str() ); | |
| 724 | ✗ | if (ierr>0) { | |
| 725 | ✗ | FEL_ERROR("error in execution of " + cmd.str() ); | |
| 726 | } | ||
| 727 | ✗ | std::stringstream filename; | |
| 728 | ✗ | filename << FelisceParam::instance().resultDir << "BVMapping.dat"; | |
| 729 | |||
| 730 | ✗ | std::ofstream outFile; | |
| 731 | ✗ | outFile.open(filename.str().c_str()); | |
| 732 | //! we write the header of the first part | ||
| 733 | ✗ | outFile<<"# Total Number of Processors"<<std::endl; | |
| 734 | ✗ | outFile<<MpiInfo::numProc()<<std::endl; | |
| 735 | ✗ | outFile<<"# Number of boundary dofs"<<std::endl; | |
| 736 | ✗ | outFile<<m_numGlobalDofInterface<<std::endl; | |
| 737 | |||
| 738 | // for each id we get the corresponding point //TODO create a function for that. | ||
| 739 | ✗ | const std::size_t iUC=m_dofPtr->getNumGlobComp(1, 0);//TODO | |
| 740 | ✗ | std::unordered_map<int,Point> id2Point; | |
| 741 | ✗ | for(auto itLabel = m_listOfBoundaryPetscDofs.at(iUC).begin(); itLabel != m_listOfBoundaryPetscDofs.at(iUC).end(); ++itLabel) { | |
| 742 | ✗ | int cpoint(0); | |
| 743 | ✗ | for(auto itDof = itLabel->second.begin(); itDof != itLabel->second.end(); ++itDof) { | |
| 744 | ✗ | id2Point[itDof->second]=m_mapLab2AllPoints.find(itLabel->first)->second[cpoint]; | |
| 745 | ✗ | cpoint++; | |
| 746 | } | ||
| 747 | } | ||
| 748 | |||
| 749 | // a second line for humans | ||
| 750 | //! We write a table like this: | ||
| 751 | //! PointsCoordinate, boundaryApplicationOrdering (BAO), boundaryPetscOrdering (BPO), volumePetscOrdering (VPO), volumeApplicationOrdering (VAO) | ||
| 752 | ✗ | outFile<<"# PointsCoordinate, boundaryApplicationOrdering (BAO), boundaryPetscOrdering (BPO), volumePetscOrdering (VPO), volumeApplicationOrdering (VAO)"<<std::endl; | |
| 753 | ✗ | for ( int cPoint(0); cPoint < m_numGlobalDofInterface; ++cPoint ) { | |
| 754 | ✗ | Point P=id2Point[m_glob2PetscVolume[cPoint]]; | |
| 755 | ✗ | felInt BPO=cPoint; AOApplicationToPetsc( m_aoInterface, 1, &BPO ); | |
| 756 | ✗ | felInt VAO=m_glob2PetscVolume[cPoint]; AOPetscToApplication( m_ao, 1, &VAO); | |
| 757 | outFile | ||
| 758 | ✗ | << P.x() <<" , " | |
| 759 | ✗ | << P.y() <<" , " | |
| 760 | ✗ | << P.z() <<" , " | |
| 761 | ✗ | << cPoint <<" , " | |
| 762 | ✗ | << BPO <<" , " | |
| 763 | ✗ | << m_glob2PetscVolume[cPoint] <<" , " | |
| 764 | ✗ | << VAO << std::endl; | |
| 765 | } | ||
| 766 | ✗ | outFile.close(); | |
| 767 | } | ||
| 768 | } | ||
| 769 | |||
| 770 | /***********************************************************************************/ | ||
| 771 | /***********************************************************************************/ | ||
| 772 | |||
| 773 | ✗ | void DofBoundary::displayBoundaryVolumeMapping() const { | |
| 774 | ✗ | if ( ! m_volumeInterfaceMappingBuilt ) | |
| 775 | ✗ | FEL_ERROR("volume interface mapping not yet build"); | |
| 776 | ✗ | const int N(25); // used only for printing information in debug | |
| 777 | typedef std::stringstream aStream; | ||
| 778 | ✗ | aStream master; | |
| 779 | ✗ | master << "m_glob2PetscVolume maps the global indeces in the application ordering (0,...," << m_numGlobalDofInterface-1 <<")" << std::endl; | |
| 780 | ✗ | master << " to the corresponding indeces of the volume in the petsc ordering (" << m_glob2PetscVolume[0]<<",...,"<<m_glob2PetscVolume[m_numGlobalDofInterface-1]<<")"<<std::endl; | |
| 781 | ✗ | PetscPrintf( MpiInfo::petscComm(), "%s",master.str().c_str() ); | |
| 782 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 783 | ✗ | if ( FelisceParam::verbose() > 3) { | |
| 784 | ✗ | aStream lineOne, lineTwo; | |
| 785 | ✗ | lineOne<<"The first "<<N<<" mapped dofs being: "<<std::endl; | |
| 786 | ✗ | for ( int ii(0); ii < N && ii < m_numGlobalDofInterface; ++ii) { | |
| 787 | ✗ | lineOne << ii << '\t' ; | |
| 788 | ✗ | lineTwo << m_glob2PetscVolume[ii] << '\t' ; | |
| 789 | } | ||
| 790 | ✗ | lineOne << std::endl; | |
| 791 | ✗ | lineTwo << std::endl; | |
| 792 | ✗ | PetscPrintf( MpiInfo::petscComm(), "%s",lineOne.str().c_str() ); | |
| 793 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 794 | ✗ | PetscPrintf( MpiInfo::petscComm(), "%s",lineTwo.str().c_str() ); | |
| 795 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 796 | } | ||
| 797 | ✗ | master.str(""); | |
| 798 | ✗ | master.clear(); | |
| 799 | ✗ | master<<"Number of local dof (for each proc) on the interface is ... "<<std::endl; | |
| 800 | ✗ | PetscPrintf(MpiInfo::petscComm(),"%s", master.str().c_str()); | |
| 801 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 802 | ✗ | std::cout<<"["<<MpiInfo::rankProc()<<"] numLoc ("<<m_numLocalDofInterface<<"/"<<m_numGlobalDofInterface<<") (local/global)"<<std::endl; | |
| 803 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 804 | ✗ | master.str(""); | |
| 805 | ✗ | master.clear(); | |
| 806 | ✗ | master<<"m_loc2GlobInterface let you go from local indexing to global ids in the application ordering"<<std::endl; | |
| 807 | ✗ | master<<"---------info are displayed only for the last proc------------"<<std::endl; | |
| 808 | ✗ | PetscPrintf( MpiInfo::petscComm(),"%s", master.str().c_str() ); | |
| 809 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 810 | ✗ | if ( MpiInfo::rankProc() == MpiInfo::numProc() - 1) { | |
| 811 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] mapping from local (0,...,"<<m_numLocalDofInterface-1<<")"<<std::endl; | |
| 812 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] The first "<< ( (N<m_numLocalDofInterface) ? N : m_numLocalDofInterface ) <<" being: "<<std::endl; | |
| 813 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] "; | |
| 814 | ✗ | for ( int ii(0); ii<N && ii<m_numLocalDofInterface; ++ii ) { | |
| 815 | ✗ | std::cout<<ii<<'\t'; | |
| 816 | } | ||
| 817 | ✗ | std::cout<<std::endl; | |
| 818 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] "; | |
| 819 | ✗ | for ( int ii(0); ii<N && ii<m_numLocalDofInterface; ++ii ) { | |
| 820 | ✗ | std::cout<<m_loc2GlobInterface[ii]<<'\t'; | |
| 821 | } | ||
| 822 | ✗ | std::cout<<std::endl; | |
| 823 | } | ||
| 824 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 825 | ✗ | master.str(""); | |
| 826 | ✗ | master.clear(); | |
| 827 | ✗ | master<<"m_loc2PetscVolume let you go from local indexing to global ids in the volume petsc ordering"<<std::endl; | |
| 828 | ✗ | master<<"---------info are displayed only for the last proc------------"<<std::endl; | |
| 829 | ✗ | PetscPrintf( MpiInfo::petscComm(),"%s", master.str().c_str() ); | |
| 830 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 831 | ✗ | if ( MpiInfo::rankProc() == MpiInfo::numProc() -1 ) { | |
| 832 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] mapping from local (0,...,"<<m_numLocalDofInterface-1<<")"<<std::endl; | |
| 833 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] The first "<<( ( N < m_numLocalDofInterface )?N:m_numLocalDofInterface)<<" being: "<<std::endl; | |
| 834 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] "; | |
| 835 | ✗ | for ( int ii(0); ii<N && ii<m_numLocalDofInterface; ++ii ) { | |
| 836 | ✗ | std::cout<<ii<<'\t'; | |
| 837 | } | ||
| 838 | ✗ | std::cout<<std::endl; | |
| 839 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] "; | |
| 840 | ✗ | for ( int ii(0); ii<N && ii<m_numLocalDofInterface; ++ii ) { | |
| 841 | ✗ | std::cout<<m_loc2PetscVolume[ii]<<'\t'; | |
| 842 | } | ||
| 843 | ✗ | std::cout<<std::endl; | |
| 844 | } | ||
| 845 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 846 | ✗ | if ( FelisceParam::verbose() > 2 ) { | |
| 847 | ✗ | master.str(""); | |
| 848 | ✗ | master.clear(); | |
| 849 | ✗ | master<<" m_aoInterface "<<std::endl; | |
| 850 | ✗ | PetscPrintf(MpiInfo::petscComm(), "%s",master.str().c_str()); | |
| 851 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 852 | ✗ | AOView(m_aoInterface,PETSC_VIEWER_STDOUT_WORLD); | |
| 853 | } | ||
| 854 | ✗ | master.str(""); | |
| 855 | ✗ | master.clear(); | |
| 856 | ✗ | master<<"m_loc2PetscInterface let you go from local indexing to global indexing petsc ordering"<<std::endl; | |
| 857 | ✗ | PetscPrintf( MpiInfo::petscComm(), "%s",master.str().c_str() ); | |
| 858 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 859 | ✗ | if ( FelisceParam::verbose() > 2 ) { | |
| 860 | ✗ | master.str(""); | |
| 861 | ✗ | master.clear(); | |
| 862 | ✗ | master<<"----------info are displayed only for the last proc"<<std::endl; | |
| 863 | ✗ | PetscPrintf( MpiInfo::petscComm(), "%s",master.str().c_str() ); | |
| 864 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 865 | ✗ | if ( MpiInfo::rankProc() == MpiInfo::numProc() - 1 ) { | |
| 866 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] mapping from local (0,...,"<<m_numLocalDofInterface-1<<")"<<std::endl; | |
| 867 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] The first "<<( ( N < m_numLocalDofInterface )?N:m_numLocalDofInterface )<<" being: "<<std::endl; | |
| 868 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] "; | |
| 869 | ✗ | for ( int ii(0); ii<N && ii<m_numLocalDofInterface; ++ii ) { | |
| 870 | ✗ | std::cout<<ii<<'\t'; | |
| 871 | } | ||
| 872 | ✗ | std::cout<<std::endl; | |
| 873 | ✗ | std::cout<<"-----------["<<MpiInfo::rankProc()<<"] "; | |
| 874 | ✗ | for ( int ii(0); ii<N && ii<m_numLocalDofInterface; ++ii ) { | |
| 875 | ✗ | std::cout<<m_loc2PetscInterface[ii]<<'\t'; | |
| 876 | } | ||
| 877 | ✗ | std::cout<<std::endl; | |
| 878 | } | ||
| 879 | ✗ | MPI_Barrier(MpiInfo::petscComm()); | |
| 880 | } | ||
| 881 | } | ||
| 882 | |||
| 883 | /***********************************************************************************/ | ||
| 884 | /***********************************************************************************/ | ||
| 885 | |||
| 886 | ✗ | void DofBoundary::unRoll( felInt iUnkComp ) { | |
| 887 | ✗ | if ( ! m_dofsHaveBeenComputed[iUnkComp] ) | |
| 888 | ✗ | FEL_ERROR("need to calculate those dofs"); | |
| 889 | ✗ | if ( m_unRolledList.size() <= std::size_t( iUnkComp ) ) { | |
| 890 | ✗ | m_unRolledList.resize(iUnkComp+1); | |
| 891 | } | ||
| 892 | ✗ | if ( m_unRolledList[iUnkComp].size() > 0 ) | |
| 893 | ✗ | return; | |
| 894 | else { | ||
| 895 | ✗ | for(auto it_label = m_listOfBoundaryPetscDofs[iUnkComp].begin(); it_label != m_listOfBoundaryPetscDofs[iUnkComp].end(); ++it_label ) { | |
| 896 | ✗ | for(auto it_point = it_label->second.begin(); it_point != it_label->second.end(); ++it_point) { | |
| 897 | ✗ | m_unRolledList[iUnkComp].insert( it_point->second ); | |
| 898 | } | ||
| 899 | } | ||
| 900 | } | ||
| 901 | } | ||
| 902 | |||
| 903 | /***********************************************************************************/ | ||
| 904 | /***********************************************************************************/ | ||
| 905 | |||
| 906 | ✗ | bool DofBoundary::isOnSurface( felInt iUnkComp, felInt petscDof ) { | |
| 907 | ✗ | if ( m_unRolledList.size() <= std::size_t( iUnkComp ) ) { | |
| 908 | ✗ | this->unRoll(iUnkComp); | |
| 909 | } | ||
| 910 | ✗ | if ( m_unRolledList[iUnkComp].size() == std::size_t(0) ) { | |
| 911 | ✗ | this->unRoll(iUnkComp); | |
| 912 | } | ||
| 913 | ✗ | auto it = m_unRolledList[iUnkComp].find( petscDof ); | |
| 914 | ✗ | return it != m_unRolledList[iUnkComp].end(); | |
| 915 | } | ||
| 916 | |||
| 917 | } | ||
| 918 |