levelSetMaskObject.cpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit is distributed in the hope that it will be useful, but WITHOUT
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *  FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
# include "levelSetMaskObject.hpp"
 
namespace bitpit {
 
LevelSetMaskObject::LevelSetMaskObject(int id, const std::unordered_set<long> &mask, const VolumeKernel &mesh )
    : LevelSetSegmentationObject(id) {
 
    std::unordered_map<long,long> meshToEnvelope ;
 
    std::unique_ptr<SurfUnstructured> segmentation = extractCellEnvelope(mask,mesh,meshToEnvelope) ;
 
    long intrIndex = meshToEnvelope.begin()->first;
    long enveIndex = meshToEnvelope.begin()->second;
 
    bool sameOrientation = sameInterfaceEnvelopeOrientation(mesh, intrIndex, *segmentation, enveIndex);
 
    auto const &interface = mesh.getInterface(intrIndex);
    long ownerId = interface.getOwner();
    bool invert = (mask.count(ownerId)==0) ;
 
    bool flip = (sameOrientation == invert);
 
    bool orientable = segmentation->adjustCellOrientation( enveIndex, flip);
    if( !orientable){
        throw std::runtime_error ("Error in LevelSetMaskObject");
    }
 
    this->setSurface(std::move(segmentation));
}
 
LevelSetMaskObject::LevelSetMaskObject(int id, const std::vector<long> &list, long interfaceId, bool invert, const VolumeKernel &mesh )
    : LevelSetSegmentationObject(id) {
 
    std::unordered_map<long,long> meshToEnvelope ;
 
    std::unique_ptr<SurfUnstructured> segmentation = extractFaceEnvelope(list,mesh,meshToEnvelope) ;
 
    long enveIndex = meshToEnvelope.at(interfaceId);
    bool sameOrientation = sameInterfaceEnvelopeOrientation(mesh, interfaceId, *segmentation, enveIndex);
 
    bool flip = (sameOrientation == invert);
 
    bool orientable = segmentation->adjustCellOrientation( enveIndex, flip);
    if( !orientable){
        throw std::runtime_error ("Error in LevelSetMaskObject");
    }
 
    this->setSurface(std::move(segmentation));
}
 
std::unique_ptr<SurfUnstructured> LevelSetMaskObject::extractCellEnvelope(const std::unordered_set<long> &mask, const VolumeKernel &mesh, std::unordered_map<long,long> &meshToEnvelope){
 
    std::vector<long> list;
 
    for( long cellIndex : mask){
        const auto &cell = mesh.getCell(cellIndex);
        const long *adjacencies = cell.getAdjacencies();
        const long *interfaceIndex = cell.getInterfaces();
 
        int interfaceCount= cell.getInterfaceCount() ;
 
        for(int i=0; i<interfaceCount; i++){
            long neigh = adjacencies[i];
 
            if(mask.count(neigh)==0){
                list.push_back(interfaceIndex[i]);
            }
        }
    }
 
    return extractFaceEnvelope(list,mesh,meshToEnvelope);
}
 
std::unique_ptr<SurfUnstructured> LevelSetMaskObject::extractFaceEnvelope(const std::vector<long> &list, const VolumeKernel &mesh, std::unordered_map<long,long> &meshToEnvelope){
 
    std::unique_ptr<SurfUnstructured> envelope;
#if BITPIT_ENABLE_MPI==1
    envelope = std::unique_ptr<SurfUnstructured>(new SurfUnstructured(mesh.getDimension()-1,mesh.getCommunicator()));
#else
    envelope = std::unique_ptr<SurfUnstructured>(new SurfUnstructured(mesh.getDimension()-1));
#endif
 
    // ====================================================================== //
    // RESIZE DATA STRUCTURES                                                 //
    // ====================================================================== //
    long nVertices(0);
    long nCells(list.size());
 
    for( long faceIndex : list){
        auto const &interface = mesh.getInterface(faceIndex);
        nVertices += interface.getVertexCount() ;
    }
 
    envelope->reserveVertices(nVertices);
    envelope->reserveCells(nCells);
 
    // ====================================================================== //
    // LOOP OVER CELLS                                                        //
    // ====================================================================== //
    std::unordered_map<long,long> vertexMap;
 
    for( long faceIndex : list){
        auto const &interface = mesh.getInterface(faceIndex);
        ConstProxyVector<long> faceVertexIds = interface.getVertexIds();
        int nFaceVertices = faceVertexIds.size();
 
        // Add face vertices to the envelope and get face
        // connectivity in the envelope
        std::unique_ptr<long[]> faceEnvelopeConnect = std::unique_ptr<long[]>(new long[nFaceVertices]);
        for (int j = 0; j < nFaceVertices; ++j) {
            long vertexId = faceVertexIds[j];
 
            // If the vertex is not yet in the envelope
            // add it.
            if (vertexMap.count(vertexId) == 0) {
                const Vertex &vertex = mesh.getVertex(vertexId);
                auto envelopeVertex = envelope->addVertex(vertex);
                vertexMap[vertexId] = envelopeVertex->getId();
            }
 
            // Update face connectivity in the envelope
            faceEnvelopeConnect[j] = vertexMap.at(vertexId);
        }
 
        // Add face to envelope
        ElementType faceType = interface.getType();
        PatchKernel::CellIterator cellItr = envelope->addCell(faceType, std::move(faceEnvelopeConnect));
        meshToEnvelope.insert({{faceIndex,cellItr.getId()}});
    }
 
    envelope->squeeze();
    envelope->initializeAdjacencies();
    envelope->initializeInterfaces();
 
    envelope->getVTK().setName("geometry_002") ;
    envelope->write() ;
 
    return envelope;
 
}
 
bool LevelSetMaskObject::sameInterfaceEnvelopeOrientation(const VolumeKernel &mesh, long faceIndex, SurfUnstructured &envelope, long enveIndex){
 
    std::array<double,3> facetNormal = envelope.evalFacetNormal(enveIndex);
    std::array<double,3> interfaceNormal = mesh.evalInterfaceNormal(faceIndex);
 
    return (dotProduct(facetNormal,interfaceNormal)>0) ;
 
}
 
}