Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

levelset: fix evaluation of levelset for surfaces that contain pixel elements #376

Merged
merged 4 commits into from
Apr 20, 2023
Merged

levelset: fix evaluation of levelset for surfaces that contain pixel elements #376

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
638ba15
patchkernel: fix orientation flip of pixel cells
andrea-iob Mar 1, 2023
b85adea
patchkernel: fix hash operator of ElementHalfItem class for pixel ele…
andrea-iob Mar 21, 2023
a1cd488
patchkernel: allow to get an anti-clockwise ordered list of facet ver…
andrea-iob Mar 1, 2023
7bdc10a
levelset: fix evaluation of segment information for segment with non-…
andrea-iob Mar 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/levelset/levelSetSegmentationObject.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,8 +189,9 @@ int SegmentationKernel::getSegmentInfo( const std::array<double,3> &pointCoords,

default:
{
ConstProxyVector<long> elementVertexIds = m_surface->getFacetOrderedVertexIds(segment);
BITPIT_CREATE_WORKSPACE(segmentVertexCoors, std::array<double BITPIT_COMMA 3>, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
m_surface->getElementVertexCoordinates(segment, segmentVertexCoors);
m_surface->getVertexCoords(elementVertexIds.size(), elementVertexIds.data(), segmentVertexCoors);
pointProjectionVector -= CGElem::projectPointPolygon( pointCoords, nSegmentVertices, segmentVertexCoors, lambda );

break;
Expand Down
6 changes: 3 additions & 3 deletions src/patchkernel/element.tpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,9 +182,9 @@ std::size_t ElementHalfItem<DerivedElement>::Hasher::operator()(const ElementHal
}
} else {
// Reversing the winding of pixel elements needs special handling. In
// pixel elements, vertices are ordered using a Z-order, wherase in
// all other elements ordering of the vertices is anti-clockwise.
if (item.m_element.getType() != ElementType::VOXEL) {
// pixel elements, vertices are ordered using the Z-order, whereas in
// all other elements vertices are ordered anti-clockwise.
if (item.m_element.getType() != ElementType::PIXEL) {
for (std::size_t i = nVertices; i > 0; --i) {
std::size_t k = (item.m_firstVertexId + i) % nVertices;
utils::hashing::hash_combine(hash, vertexIds[k]);
Expand Down
216 changes: 136 additions & 80 deletions src/patchkernel/surface_kernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,9 +279,9 @@ double SurfaceKernel::evalCellArea(long id) const
double coeff = 0.25;
double area = 0.0;
for (int i = 0; i < nvert; ++i) {
int prevVertex = getOrderedLocalVertexIds(*cell_, (nvert + i - 1) % nvert);
int vertex = getOrderedLocalVertexIds(*cell_, i);
int nextVertex = getOrderedLocalVertexIds(*cell_, (i + 1) % nvert);
int prevVertex = getFacetOrderedLocalVertex(*cell_, (nvert + i - 1) % nvert);
int vertex = getFacetOrderedLocalVertex(*cell_, i);
int nextVertex = getFacetOrderedLocalVertex(*cell_, (i + 1) % nvert);

const std::array<double, 3> &prevVertexCoords = m_vertices[cellVertexIds[prevVertex]].getCoords();
const std::array<double, 3> &vertexCoords = m_vertices[cellVertexIds[vertex]].getCoords();
Expand Down Expand Up @@ -457,8 +457,8 @@ double SurfaceKernel::evalAngleAtVertex(long id, int vertex) const
ConstProxyVector<long> cellVertexIds = cell.getVertexIds();
int nCellVertices = cellVertexIds.size();

int prevVertex = getOrderedLocalVertexIds(cell, (vertex - 1 + nCellVertices) % nCellVertices);
int nextVertex = getOrderedLocalVertexIds(cell, (vertex + 1) % nCellVertices);
int prevVertex = getFacetOrderedLocalVertex(cell, (vertex - 1 + nCellVertices) % nCellVertices);
int nextVertex = getFacetOrderedLocalVertex(cell, (vertex + 1) % nCellVertices);

const std::array<double, 3> &prevVertexCoords = m_vertices[cellVertexIds[prevVertex]].getCoords();
const std::array<double, 3> &vertexCoords = m_vertices[cellVertexIds[vertex]].getCoords();
Expand Down Expand Up @@ -651,9 +651,9 @@ std::array<double, 3> SurfaceKernel::evalFacetNormal(long id, const std::array<d
int i, nvert = cellVertexIds.size();
double coeff = 1.0/double(nvert);
for (i = 0; i < nvert; ++i) {
int prevVertex = getOrderedLocalVertexIds(*cell_, (nvert + i - 1) % nvert);
int vertex = getOrderedLocalVertexIds(*cell_, i);
int nextVertex = getOrderedLocalVertexIds(*cell_, (i + 1) % nvert);
int prevVertex = getFacetOrderedLocalVertex(*cell_, (nvert + i - 1) % nvert);
int vertex = getFacetOrderedLocalVertex(*cell_, i);
int nextVertex = getFacetOrderedLocalVertex(*cell_, (i + 1) % nvert);

const std::array<double, 3> &prevVertexCoords = m_vertices[cellVertexIds[prevVertex]].getCoords();
const std::array<double, 3> &vertexCoords = m_vertices[cellVertexIds[vertex]].getCoords();
Expand Down Expand Up @@ -1158,16 +1158,16 @@ bool SurfaceKernel::haveSameOrientation(const Cell &cell_A, int face_A, const Ce
// orientation, if the vertices appear in reversed order, the facets have
// the same orientation.
for (std::size_t i = 0; i < nCellVertices_A; ++i) {
long vertexId_A = cellVertexIds_A[getOrderedLocalVertexIds(cell_A, i)];
long vertexId_A = cellVertexIds_A[getFacetOrderedLocalVertex(cell_A, i)];
for (std::size_t j = 0; j < nCellVertices_B; ++j) {
long vertexId_B = cellVertexIds_B[getOrderedLocalVertexIds(cell_B, j)];
long vertexId_B = cellVertexIds_B[getFacetOrderedLocalVertex(cell_B, j)];
if (vertexId_A == vertexId_B) {
if (cellDimension == 2) {
long previousVertexId_A = cellVertexIds_A[getOrderedLocalVertexIds(cell_A, (i - 1 + nCellVertices_A) % nCellVertices_A)];
long previousVertexId_B = cellVertexIds_B[getOrderedLocalVertexIds(cell_B, (j - 1 + nCellVertices_B) % nCellVertices_B)];
long previousVertexId_A = cellVertexIds_A[getFacetOrderedLocalVertex(cell_A, (i - 1 + nCellVertices_A) % nCellVertices_A)];
long previousVertexId_B = cellVertexIds_B[getFacetOrderedLocalVertex(cell_B, (j - 1 + nCellVertices_B) % nCellVertices_B)];

long nextVertexId_A = cellVertexIds_A[getOrderedLocalVertexIds(cell_A, (i + 1) % nCellVertices_A)];
long nextVertexId_B = cellVertexIds_B[getOrderedLocalVertexIds(cell_B, (j + 1) % nCellVertices_B)];
long nextVertexId_A = cellVertexIds_A[getFacetOrderedLocalVertex(cell_A, (i + 1) % nCellVertices_A)];
long nextVertexId_B = cellVertexIds_B[getFacetOrderedLocalVertex(cell_B, (j + 1) % nCellVertices_B)];

if (nextVertexId_A == nextVertexId_B || previousVertexId_A == previousVertexId_B) {
return false;
Expand Down Expand Up @@ -1318,7 +1318,8 @@ void SurfaceKernel::flipCellOrientation(long id)
// Flipped faces 3 2 1 0 4
//
// We have reversed the order of the vertices, this means that we have
// to invert the order of all the faces but the lust one.
// to invert the order of all the faces but the lust one. Pixel elements
// should be handled separately.
//
// For one-dimensional elements faces are defined on the vertices itself:
//
Expand All @@ -1330,75 +1331,85 @@ void SurfaceKernel::flipCellOrientation(long id)
//
// We have reversed the order of the vertices, this means that we have
// to invert the order of all the faces.
int nCellFaces = cell.getFaceCount();
int nCellAdjacencies = cell.getAdjacencyCount();
int nCellInterfaces = cell.getInterfaceCount();

int flipOffset;
if (cell.getDimension() == 1) {
flipOffset = 0;
} else {
flipOffset = 1;
}
if (cell.getType() != ElementType::PIXEL) {
int nCellFaces = cell.getFaceCount();
int nCellAdjacencies = cell.getAdjacencyCount();
int nCellInterfaces = cell.getInterfaceCount();

FlatVector2D<long> flippedCellAdjacencies;
FlatVector2D<long> flippedCellInterfaces;
flippedCellAdjacencies.reserve(nCellFaces, nCellAdjacencies);
flippedCellInterfaces.reserve(nCellFaces, nCellInterfaces);
int flipOffset;
if (cell.getDimension() == 1) {
flipOffset = 0;
} else {
flipOffset = 1;
}

for (int i = 0; i < nCellFaces - flipOffset; ++i) {
int nFaceAdjacencies = cell.getAdjacencyCount(nCellFaces - i - 1 - flipOffset);
const long *faceAdjacencies = cell.getAdjacencies(nCellFaces - i - 1 - flipOffset);
flippedCellAdjacencies.pushBack(nFaceAdjacencies, faceAdjacencies);
FlatVector2D<long> flippedCellAdjacencies;
FlatVector2D<long> flippedCellInterfaces;
flippedCellAdjacencies.reserve(nCellFaces, nCellAdjacencies);
flippedCellInterfaces.reserve(nCellFaces, nCellInterfaces);

int nFaceInterfaces = cell.getInterfaceCount(nCellFaces - i - 1 - flipOffset);
const long *faceInterfaces = cell.getInterfaces(nCellFaces - i - 1 - flipOffset);
flippedCellInterfaces.pushBack(nFaceInterfaces, faceInterfaces);
}
for (int i = 0; i < nCellFaces - flipOffset; ++i) {
int nFaceAdjacencies = cell.getAdjacencyCount(nCellFaces - i - 1 - flipOffset);
const long *faceAdjacencies = cell.getAdjacencies(nCellFaces - i - 1 - flipOffset);
flippedCellAdjacencies.pushBack(nFaceAdjacencies, faceAdjacencies);

if (flipOffset == 1) {
int nLastFaceAdjacencies = cell.getAdjacencyCount(nCellFaces - 1);
const long *lastFaceAdjacencies = cell.getAdjacencies(nCellFaces - 1);
flippedCellAdjacencies.pushBack(nLastFaceAdjacencies, lastFaceAdjacencies);
int nFaceInterfaces = cell.getInterfaceCount(nCellFaces - i - 1 - flipOffset);
const long *faceInterfaces = cell.getInterfaces(nCellFaces - i - 1 - flipOffset);
flippedCellInterfaces.pushBack(nFaceInterfaces, faceInterfaces);
}

int nLastFaceInterfaces = cell.getInterfaceCount(nCellFaces - 1);
const long *lastFaceInterfaces = cell.getInterfaces(nCellFaces - 1);
flippedCellInterfaces.pushBack(nLastFaceInterfaces, lastFaceInterfaces);
}
if (flipOffset == 1) {
int nLastFaceAdjacencies = cell.getAdjacencyCount(nCellFaces - 1);
const long *lastFaceAdjacencies = cell.getAdjacencies(nCellFaces - 1);
flippedCellAdjacencies.pushBack(nLastFaceAdjacencies, lastFaceAdjacencies);

// Set flipped adjacencies
if (nCellAdjacencies > 0) {
cell.setAdjacencies(std::move(flippedCellAdjacencies));
}
int nLastFaceInterfaces = cell.getInterfaceCount(nCellFaces - 1);
const long *lastFaceInterfaces = cell.getInterfaces(nCellFaces - 1);
flippedCellInterfaces.pushBack(nLastFaceInterfaces, lastFaceInterfaces);
}

// Set flipped interfaces
//
// After setting the flipped interfaces, we also need to update the face
// associated with the interfaces.
if (nCellInterfaces > 0) {
cell.setInterfaces(std::move(flippedCellInterfaces));
// Set flipped adjacencies
if (nCellAdjacencies > 0) {
cell.setAdjacencies(std::move(flippedCellAdjacencies));
}

const long *cellInterfaces = cell.getInterfaces();
int nCellInterfaces = cell.getInterfaceCount();
for (int i = 0; i < nCellInterfaces; ++i) {
long interfaceId = cellInterfaces[i];
Interface &interface = getInterface(interfaceId);

long owner = interface.getOwner();
if (owner == id) {
int ownerFace = interface.getOwnerFace();
if (ownerFace != nCellFaces - 1) {
ownerFace = nCellFaces - 2 - ownerFace;
interface.setOwner(id, ownerFace);
}
} else {
int neighFace = interface.getNeighFace();
if (neighFace != nCellFaces - 1) {
neighFace = nCellFaces - 2 - neighFace;
interface.setNeigh(id,neighFace);
// Set flipped interfaces
//
// After setting the flipped interfaces, we also need to update the face
// associated with the interfaces.
if (nCellInterfaces > 0) {
cell.setInterfaces(std::move(flippedCellInterfaces));

const long *cellInterfaces = cell.getInterfaces();
int nCellInterfaces = cell.getInterfaceCount();
for (int i = 0; i < nCellInterfaces; ++i) {
long interfaceId = cellInterfaces[i];
Interface &interface = getInterface(interfaceId);

long owner = interface.getOwner();
if (owner == id) {
int ownerFace = interface.getOwnerFace();
if (ownerFace != nCellFaces - 1) {
ownerFace = nCellFaces - 2 - ownerFace;
interface.setOwner(id, ownerFace);
}
} else {
int neighFace = interface.getNeighFace();
if (neighFace != nCellFaces - 1) {
neighFace = nCellFaces - 2 - neighFace;
interface.setNeigh(id,neighFace);
}
}
}
}
} else {
long *adjacencies = cell.getAdjacencies();
std::swap(adjacencies[0], adjacencies[1]);
std::swap(adjacencies[2], adjacencies[3]);

long *interfaces = cell.getInterfaces();
std::swap(interfaces[0], interfaces[1]);
std::swap(interfaces[2], interfaces[3]);
}
}

Expand Down Expand Up @@ -1651,17 +1662,61 @@ bool SurfaceKernel::compareSelectedTypes(unsigned short mask_, ElementType type_
}

/*!
* Get the local index of the n-th vertex in the anti-clockwise ordered list of
* vertex ids.
Get the anti-clockwise ordered list of vertex for the specified facet.

\param facet is the facet
\result The anti-clockwise ordered list of vertex for the specified facet.
*/
ConstProxyVector<long> SurfaceKernel::getFacetOrderedVertexIds(const Cell &facet) const
{
ConstProxyVector<long> vertexIds = facet.getVertexIds();
if (areFacetVerticesOrdered(facet)) {
return vertexIds;
}

std::size_t nVertices = vertexIds.size();
ConstProxyVector<long> orderedVertexIds(ConstProxyVector<long>::INTERNAL_STORAGE, nVertices);
ConstProxyVector<long>::storage_pointer orderedVertexIdsStorage = orderedVertexIds.storedData();
for (std::size_t k = 0; k < nVertices; ++k) {
int vertex = getFacetOrderedLocalVertex(facet, k);
orderedVertexIdsStorage[k] = vertexIds[vertex];
}

return orderedVertexIds;
}

/*!
* Check if the vertices of the specified facet are anti-clockwise ordered.
*
* \param[in] cell is the cell
* \param[in] n is the index of the requested vertex
* \result The the local index of the n-th vertex in the anti-clockwise ordered
* \param[in] facet is the facet
* \result Return true if the vertices of the specified facet are anti-clockwise ordered,
* false otherwise.
*/
bool SurfaceKernel::areFacetVerticesOrdered(const Cell &facet) const
{
switch (facet.getType()) {

case ElementType::PIXEL:
return false;

default:
return true;

}
}

/*!
* Get the local index of the vertex occupying the n-th position in the anti-clockwise ordered
* list of vertex ids.
*
* \param[in] facet is the facet
* \param[in] n is the requested position
* \result The local index of the vertex occupying the n-th position in the anti-clockwise
* ordered list of vertex ids.
*/
int SurfaceKernel::getOrderedLocalVertexIds(const Cell &cell, long n) const
int SurfaceKernel::getFacetOrderedLocalVertex(const Cell &facet, std::size_t n) const
{
switch (cell.getType()) {
switch (facet.getType()) {

case ElementType::PIXEL:
if (n == 2) {
Expand All @@ -1673,6 +1728,7 @@ int SurfaceKernel::getOrderedLocalVertexIds(const Cell &cell, long n) const
}

default:
assert(areFacetVerticesOrdered(facet));
return n;

}
Expand Down
5 changes: 4 additions & 1 deletion src/patchkernel/surface_kernel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,8 @@ class SurfaceKernel : public PatchKernel {
bool adjustCellOrientation(long id, bool invert = false);
void flipCellOrientation(long id);

ConstProxyVector<long> getFacetOrderedVertexIds(const Cell &facet) const;

void displayQualityStats(std::ostream&, unsigned int padding = 0) const;
std::vector<double> computeHistogram(eval_f_ funct_, std::vector<double> &bins, long &count, int n_intervals = 8, unsigned short mask = SELECT_ALL) const;

Expand All @@ -97,7 +99,8 @@ class SurfaceKernel : public PatchKernel {
SurfaceKernel(int id, int dimension, bool expert);
#endif

int getOrderedLocalVertexIds(const Cell &cell, long n) const;
bool areFacetVerticesOrdered(const Cell &facet) const;
int getFacetOrderedLocalVertex(const Cell &facet, std::size_t n) const;

};

Expand Down