Loading...
Searching...
No Matches
levelSetSegmentationObject.cpp
1/*---------------------------------------------------------------------------*\
2 *
3 * bitpit
4 *
5 * Copyright (C) 2015-2021 OPTIMAD engineering Srl
6 *
7 * -------------------------------------------------------------------------
8 * License
9 * This file is part of bitpit.
10 *
11 * bitpit is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU Lesser General Public License v3 (LGPL)
13 * as published by the Free Software Foundation.
14 *
15 * bitpit is distributed in the hope that it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 * License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public License
21 * along with bitpit. If not, see <http://www.gnu.org/licenses/>.
22 *
23\*---------------------------------------------------------------------------*/
24
25#include "levelSetSegmentationObject.hpp"
26
27#include "bitpit_CG.hpp"
28#include "levelSetObject.hpp"
29
30namespace bitpit {
31
43const double LevelSetSegmentationSurfaceInfo::DEFAULT_FEATURE_ANGLE = 2. * BITPIT_PI;
44
48LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo()
49 : m_surface(nullptr),
50 m_featureAngle(0),
51 m_surfaceSmoothing(LevelSetSurfaceSmoothing::LOW_ORDER)
52{
53}
54
58LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(const LevelSetSegmentationSurfaceInfo &other)
59 : m_surface(other.m_surface),
60 m_featureAngle(other.m_featureAngle),
61 m_surfaceSmoothing(other.m_surfaceSmoothing),
62 m_segmentVertexOffset(other.m_segmentVertexOffset),
63 m_segmentNormalsValid(other.m_segmentNormalsValid),
64 m_segmentNormalsStorage(other.m_segmentNormalsStorage),
65 m_unlimitedVertexNormalsValid(other.m_unlimitedVertexNormalsValid),
66 m_unlimitedVertexNormalsStorage(other.m_unlimitedVertexNormalsStorage),
67 m_limitedSegmentVertexNormalValid(other.m_limitedSegmentVertexNormalValid),
68 m_limitedSegmentVertexNormalStorage(other.m_limitedSegmentVertexNormalStorage)
69{
70 if (other.m_ownedSurface) {
71 m_ownedSurface = std::unique_ptr<SurfUnstructured>(new SurfUnstructured(*(other.m_ownedSurface)));
72 } else {
73 m_ownedSurface = nullptr;
74 }
75
76 if (other.m_searchTree) {
77 m_searchTree = std::unique_ptr<SurfaceSkdTree>(new SurfaceSkdTree(m_surface));
78 m_searchTree->build();
79 } else {
80 m_searchTree = nullptr;
81 }
82}
83
89LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(LevelSetSurfaceSmoothing surfaceSmoothing)
90 : m_surface(nullptr),
91 m_featureAngle(0),
92 m_surfaceSmoothing(surfaceSmoothing)
93{
94}
95
103LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing)
104 : m_surfaceSmoothing(surfaceSmoothing)
105{
106 setSurface(std::move(surface), featureAngle);
107}
108
116LevelSetSegmentationSurfaceInfo::LevelSetSegmentationSurfaceInfo(const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing)
117{
118 setSurface(surface, featureAngle, surfaceSmoothing);
119}
120
125const SurfUnstructured & LevelSetSegmentationSurfaceInfo::getSurface() const {
126 return *m_surface;
127}
128
134void LevelSetSegmentationSurfaceInfo::setSurface(std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle){
135
136 m_ownedSurface = std::move(surface);
137 setSurface(m_ownedSurface.get(), featureAngle);
138}
139
146void LevelSetSegmentationSurfaceInfo::setSurface(const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing){
147
148 // Check if adjacencies are built
149 if (surface->getAdjacenciesBuildStrategy() == SurfUnstructured::ADJACENCIES_NONE) {
150 throw std::runtime_error ("Segmentation needs adjacencies!") ;
151 }
152
153 // Surface information
154 m_surface = surface;
155 m_featureAngle = featureAngle;
156 m_surfaceSmoothing = surfaceSmoothing;
157
158 // Segment vertices information
159 m_segmentVertexOffset.unsetKernel();
160 m_segmentVertexOffset.setStaticKernel(&m_surface->getCells());
161
162 std::size_t nTotalSegmentVertices = 0;
163 for (auto itr = m_segmentVertexOffset.begin(); itr != m_segmentVertexOffset.end(); ++itr) {
164 *itr = nTotalSegmentVertices;
165 nTotalSegmentVertices += m_surface->getCells().rawAt(itr.getRawIndex()).getVertexCount();
166 }
167
168 // Normals
169 m_segmentNormalsValid.unsetKernel();
170 m_segmentNormalsValid.setStaticKernel(&m_surface->getCells());
171 m_segmentNormalsValid.fill(false);
172 m_segmentNormalsStorage.unsetKernel();
173 m_segmentNormalsStorage.setStaticKernel(&m_surface->getCells());
174
175 m_unlimitedVertexNormalsValid.unsetKernel();
176 m_unlimitedVertexNormalsValid.setStaticKernel(&m_surface->getVertices());
177 m_unlimitedVertexNormalsValid.fill(false);
178 m_unlimitedVertexNormalsStorage.unsetKernel();
179 m_unlimitedVertexNormalsStorage.setStaticKernel(&m_surface->getVertices());
180
181 m_limitedSegmentVertexNormalValid.assign(nTotalSegmentVertices, false);
182
183 // Initialize search tree
184 m_searchTree = std::unique_ptr<SurfaceSkdTree>(new SurfaceSkdTree(surface));
185 m_searchTree->build();
186}
187
193void LevelSetSegmentationSurfaceInfo::setSurface(const SurfUnstructured *surface, double featureAngle){
194 setSurface(surface, featureAngle, LevelSetSurfaceSmoothing::LOW_ORDER);
195}
196
201const SurfaceSkdTree & LevelSetSegmentationSurfaceInfo::getSearchTree() const {
202 return *m_searchTree;
203}
204
209double LevelSetSegmentationSurfaceInfo::getFeatureAngle() const {
210 return m_featureAngle;
211}
212
221
231double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array<double, 3> &point,
232 const SegmentConstIterator &segmentItr,
233 bool signedDistance,
234 std::array<double, 3> *distanceVector) const
235{
236 // Project the point on the surface and evaluate the point-projection vector
237 int nSegmentVertices = segmentItr->getVertexCount();
238 BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
239 std::array<double, 3> pointProjection = evalProjection(point, segmentItr, lambda);
240 (*distanceVector) = point - pointProjection;
241
242 // Evaluate unsigned distance
243 double unsignedDistance = norm2(*distanceVector);
244 if (!signedDistance) {
245 return unsignedDistance;
246 }
247
248 // Signed distance
249 //
250 // If the sign is null and the point doesn't lie on the segmentation, it lies on the normal
251 // plane. This case is not supported, because it would require to evaluate the sign taking
252 // into account the the curvature of the surface.
253 std::array<double, 3> pseudoNormal = computePseudoNormal(segmentItr, lambda);
254 double pointProjectionNormalComponent = dotProduct(*distanceVector, pseudoNormal);
255
256 double distanceTolerance = m_surface->getTol();
257 if (utils::DoubleFloatingEqual()(pointProjectionNormalComponent, 0., distanceTolerance, distanceTolerance)) {
258 bool pointOnSegmentation = utils::DoubleFloatingEqual()(unsignedDistance, 0., distanceTolerance, distanceTolerance);
259 if (!pointOnSegmentation) {
260 throw std::runtime_error("Unable to evaluate point sign: the point lies on the normal plane!");
261 }
262 }
263
264 return sign(pointProjectionNormalComponent) * unsignedDistance;
265}
266
275double LevelSetSegmentationSurfaceInfo::evalDistance(const std::array<double, 3> &point,
276 const SegmentConstIterator &segmentItr,
277 bool signedDistance) const
278{
279 std::array<double, 3> distanceVector;
280 return evalDistance(point, segmentItr, signedDistance, &distanceVector);
281}
282
290std::array<double, 3> LevelSetSegmentationSurfaceInfo::evalDistanceVector(const std::array<double, 3> &point,
291 const SegmentConstIterator &segmentItr) const
292{
293 int nSegmentVertices = segmentItr->getVertexCount();
294 BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
295 std::array<double, 3> pointProjection = evalProjection(point, segmentItr, lambda);
296
297 return (point - pointProjection);
298}
299
307std::array<double, 3> LevelSetSegmentationSurfaceInfo::evalNormal(const std::array<double, 3> &point,
308 const SegmentConstIterator &segmentItr) const
309{
310 // Project the point on the surface and evaluate the point-projection vector
311 std::array<double, 3> projectionPoint;
312 std::array<double, 3> projectionNormal;
313 evalProjection(point, segmentItr, &projectionPoint, &projectionNormal);
314
315 return projectionNormal;
316}
317
325std::array<double,3> LevelSetSegmentationSurfaceInfo::evalPseudoNormal(const SegmentConstIterator &segmentItr,
326 const double *lambda ) const
327{
328 return computePseudoNormal(segmentItr, lambda);
329}
330
345void LevelSetSegmentationSurfaceInfo::evalProjectionOnVertex(const std::array<double,3> &point,
346 const SegmentConstIterator &segmentItr,
347 std::array<double, 3> *projectionPoint,
348 std::array<double, 3> *projectionNormal) const
349{
350 BITPIT_UNUSED(point);
351
352 // Get segment
353 const Cell &segment = *segmentItr;
354 assert(segment.getType() == ElementType::VERTEX);
355
356 // Get vertex id
357 ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
358 long id = segmentVertexIds[0];
359
360 // Compute projection point and normal
361 (*projectionPoint) = m_surface->getVertexCoords(id);
362 (*projectionNormal) = {0., 0., 0.};
363}
364
415void LevelSetSegmentationSurfaceInfo::evalHighOrderProjectionOnLine(const std::array<double,3> &point,
416 const SegmentConstIterator &segmentItr,
417 std::array<double, 3> *projectionPoint,
418 std::array<double, 3> *projectionNormal) const
419{
420 // Get segment
421 const Cell &segment = *segmentItr;
422 assert(segment.getType() == ElementType::LINE);
423
424 // Get segment vertices
425 ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
426
427 const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
428 const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
429
430 // Get projection point on segment
431 std::array<double,2> t;
432 std::array<double, 3> point_s = CGElem::projectPointSegment(point, point0, point1, &(t[0]));
433
434 // Early return if point_s consides with segment's node
435 double distanceTolerance = m_surface->getTol();
436 for (int i = 0; i < 2; ++i) {
437 if (utils::DoubleFloatingEqual()(t[i], 1., distanceTolerance, distanceTolerance)) {
438 (*projectionPoint) = m_surface->getVertexCoords(segmentVertexIds[i]);
439 (*projectionNormal) = computeSegmentVertexNormal(segmentItr, i, true);
440 return;
441 }
442 }
443
444 // Get normal on segment
445 std::array<double, 3> facetNormal = computeSegmentNormal(segmentItr);
446 std::array<double, 3> normal_s = point - point_s;
447
448 double distance = norm2(normal_s);
449 if (utils::DoubleFloatingEqual()(distance, 0., distanceTolerance, distanceTolerance)) {
450 normal_s = facetNormal;
451 } else {
452 normal_s /= distance;
453 if (dotProduct(normal_s, facetNormal) < 0.0) {
454 normal_s *= -1.0;
455 }
456 }
457
458 // Get normals on segment's vertices
459 std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, 0, true);
460 std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, 1, true);
461
462 // Compute lambda coefficients
463 std::array<double, 3> edge = point1 - point0;
464 double lambda01 = -dotProduct(edge, normal0) / dotProduct(normal_s, normal0);
465 double lambda10 = dotProduct(edge, normal1) / dotProduct(normal_s, normal1);
466
467 // Eval polynomial parameterizing the curve
468 double product1 = t[0] * t[1];
469 double product2 = lambda01 * t[0] + lambda10 * t[1];
470 double f = product1 * product2;
471 double df = product1 * (lambda01 - lambda10) + product2 * (t[1] - t[0]);
472
473 // Eval projection point on curve
474 (*projectionPoint) = point_s + f * normal_s;
475
476 // Eval normal vector on curve
477 double length = norm2(edge);
478 (*projectionNormal) = df * edge + length * length * normal_s;
479 (*projectionNormal) /= norm2((*projectionNormal));
480
481 if (dotProduct(normal_s, (*projectionNormal)) < 0.0) {
482 (*projectionNormal) *= -1.0;
483 }
484}
485
616void LevelSetSegmentationSurfaceInfo::evalHighOrderProjectionOnTriangle(const std::array<double,3> &point,
617 const SegmentConstIterator &segmentItr,
618 std::array<double, 3> *projectionPoint,
619 std::array<double, 3> *projectionNormal) const
620{
621 // Get segment
622 const Cell &segment = *segmentItr;
623 assert(segment.getType() == ElementType::TRIANGLE);
624
625 // Get segment vertices
626 ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
627
628 const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
629 const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
630 const std::array<double,3> &point2 = m_surface->getVertexCoords(segmentVertexIds[2]);
631
632 // Get projection point on segment
633 std::array<double,3> tau;
634 std::array<double, 3> point_s = CGElem::projectPointTriangle(point, point0, point1, point2, &(tau[0]));
635 std::array<double, 3> d0_point_s = point0 - point2;
636 std::array<double, 3> d1_point_s = point1 - point2;
637
638 // Early return if point_s consides with segment's node
639 double distanceTolerance = m_surface->getTol();
640 for (int i = 0; i < 3; ++i) {
641 if (utils::DoubleFloatingEqual()(tau[i], 1., distanceTolerance, distanceTolerance)) {
642 (*projectionPoint) = m_surface->getVertexCoords(segmentVertexIds[i]);
643 (*projectionNormal) = computeSegmentVertexNormal(segmentItr, i, true);
644 return;
645 }
646 }
647
648 // Get normal on segment
649 std::array<double, 3> facetNormal = computeSegmentNormal(segmentItr);
650 std::array<double, 3> normal_s = point - point_s;
651
652 double distance = norm2(normal_s);
653 if (utils::DoubleFloatingEqual()(distance, 0., distanceTolerance, distanceTolerance)) {
654 normal_s = facetNormal;
655 } else {
656 normal_s /= distance;
657 if (dotProduct(normal_s, facetNormal) < 0.0) {
658 normal_s *= -1.0;
659 }
660 }
661
662 // Compute projection point and normal at each edge "e"
663 std::array<std::array<double, 3>,3 > edgePoint_s;
664 std::array< std::array<double, 3>, 3 > d0_edgePoint_s;
665 std::array< std::array<double, 3>, 3 > d1_edgePoint_s;
666 std::array< std::array<double, 3>, 3 > edgePoint_p;
667 std::array< std::array<double, 3>, 3 > d0_edgePoint_p;
668 std::array< std::array<double, 3>, 3 > d1_edgePoint_p;
669 std::array< std::array<double, 3>, 3 > edgeNormal_p;
670 std::array< std::array<double, 3>, 3 > d0_edgeNormal_p;
671 std::array< std::array<double, 3>, 3 > d1_edgeNormal_p;
672 for (int e = 0; e < 3; ++e) {
673 // Get segment
674 const Cell &segment = *segmentItr;
675
676 // Get edge's local and global node ids
677 ConstProxyVector<int> edgeLocalVertexIds = segment.getFaceLocalVertexIds(e);
678 int localId0 = edgeLocalVertexIds[0];
679 int localId1 = edgeLocalVertexIds[1];
680 int id0 = segment.getVertexId(localId0);
681 int id1 = segment.getVertexId(localId1);
682
683 // Get vertex information
684 const std::array<double,3> &node0 = m_surface->getVertexCoords(id0);
685 const std::array<double,3> &node1 = m_surface->getVertexCoords(id1);
686
687 std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, localId0, true);
688 std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, localId1, true);
689
690 // Get edge normal
691 std::array<double, 3> edgeNormal_s = computeSegmentEdgeNormal(segmentItr, e, true);
692
693 // Compute projection coefficients
694 std::array<double, 3> edge = node1 - node0;
695 double p01 = dotProduct(point_s - node0, edge) / dotProduct(edge, edge);
696 double d0_p01;
697 double d1_p01;
698 if (p01 > 1.0) {
699 p01 = 1.0;
700 d0_p01 = 0.0;
701 d1_p01 = 0.0;
702 } else if (p01 < 0.0) {
703 p01 = 0.0;
704 d0_p01 = 0.0;
705 d1_p01 = 0.0;
706 } else {
707 d0_p01 = dotProduct(d0_point_s - node0, edge) / dotProduct(edge, edge);
708 d1_p01 = dotProduct(d1_point_s - node0, edge) / dotProduct(edge, edge);
709 }
710 double p10 = 1.0 - p01;
711 double d0_p10 = - d0_p01;
712 double d1_p10 = - d1_p01;
713
714 // Compute lambda coefficients
715 double lambda01 = -dotProduct(edge, normal0) / dotProduct(edgeNormal_s, normal0);
716 double lambda10 = dotProduct(edge, normal1) / dotProduct(edgeNormal_s, normal1);
717
718 // Edge projection point
719 double product1 = p10 * p01;
720 double d_product1 = p01 - p10; // derivative wrt p10
721
722 double product2 = lambda01 * p10 + lambda10 * p01;
723 double d_product2 = lambda01 - lambda10;
724
725 double f01 = product1 * product2;
726 double d_f01 = d_product1 * product2 + product1 * d_product2;
727 double d_d_f01 = 2.0 * (d_product1 * d_product2 - product2);
728
729 edgePoint_s[e] = p10 * node0 + p01 * node1;
730 std::array<double, 3> d_edgePoint_s = node0 - node1;
731
732 d0_edgePoint_s[e] = d_edgePoint_s * d0_p10;
733 d1_edgePoint_s[e] = d_edgePoint_s * d1_p10;
734
735 edgePoint_p[e] = edgePoint_s[e] + f01 * edgeNormal_s;
736 std::array<double, 3> d_edgePoint_p = d_edgePoint_s + d_f01 * edgeNormal_s;
737
738 d0_edgePoint_p[e] = d_edgePoint_p * d0_p10;
739 d1_edgePoint_p[e] = d_edgePoint_p * d1_p10;
740
741 // Edge surface normal
742 double length = norm2(edge);
743 edgeNormal_p[e] = d_f01 * edge + length * length * edgeNormal_s;
744 std::array<double, 3> d_edgeNormal_p = d_d_f01 * edge;
745
746 double norm = norm2(edgeNormal_p[e]);
747 double d_norm = dotProduct(edgeNormal_p[e], d_edgeNormal_p) / norm;
748 edgeNormal_p[e] /= norm;
749 d_edgeNormal_p = (d_edgeNormal_p - edgeNormal_p[e] * d_norm) / norm;
750
751 d0_edgeNormal_p[e] = d_edgeNormal_p * d0_p10;
752 d1_edgeNormal_p[e] = d_edgeNormal_p * d1_p10;
753 }
754
755 // Project surface point to the intermediator polygon
756 double sum_t = 0.0;
757 double d0_sum_t = 0.0;
758 double d1_sum_t = 0.0;
759 std::array<double, 3> t;
760 std::array<double, 3> d0_t;
761 std::array<double, 3> d1_t;
762 for (int i = 0; i < 3; i ++) { // loop in vertices
763 int j = (i + 1) % 3;
764 int k = (i + 2) % 3;
765 std::array<double, 3> v1 = edgePoint_s[j] - point_s;
766 std::array<double, 3> d0_v1 = d0_edgePoint_s[j]- d0_point_s;
767 std::array<double, 3> d1_v1 = d1_edgePoint_s[j]- d1_point_s;
768
769 std::array<double, 3> v2 = edgePoint_s[k] - point_s;
770 std::array<double, 3> d0_v2 = d0_edgePoint_s[k]- d0_point_s;
771 std::array<double, 3> d1_v2 = d1_edgePoint_s[k]- d1_point_s;
772
773 std::array<double, 3> normal = crossProduct(v1, v2);
774 std::array<double, 3> d0_normal = crossProduct(d0_v1, v2) + crossProduct(v1, d0_v2);
775 std::array<double, 3> d1_normal = crossProduct(d1_v1, v2) + crossProduct(v1, d1_v2);
776
777 t[i] = norm2(normal);
778 d0_t[i] = dotProduct(normal, d0_normal) / t[i];
779 d1_t[i] = dotProduct(normal, d1_normal) / t[i];
780 sum_t += t[i];
781 d0_sum_t += d0_t[i];
782 d1_sum_t += d1_t[i];
783 }
784
785 std::array<double, 3> point_m = {0.0, 0.0, 0.0};
786 std::array<double, 3> d0_point_m = {0.0, 0.0, 0.0};
787 std::array<double, 3> d1_point_m = {0.0, 0.0, 0.0};
788 for (int i = 0; i < 3; i ++) {
789 t[i] /= sum_t;
790 d0_t[i] = (d0_t[i] - t[i] * d0_sum_t) / sum_t;
791 d1_t[i] = (d1_t[i] - t[i] * d1_sum_t) / sum_t;
792 point_m += t[i] * edgePoint_p[i];
793 d0_point_m += d0_t[i] * edgePoint_p[i] + t[i] * d0_edgePoint_p[i];
794 d1_point_m += d1_t[i] * edgePoint_p[i] + t[i] * d1_edgePoint_p[i];
795 }
796
797 // Early return if point_s coinsides with a triangle
798 // edge (or a node of the intermediate triangle)
799 for (int i = 0; i < 3; ++i) {
800 if (utils::DoubleFloatingEqual()(t[i], 1., distanceTolerance, distanceTolerance)) {
801 (*projectionPoint) = edgePoint_p[i];
802 (*projectionNormal) = edgeNormal_p[i];
803 return;
804 }
805 }
806
807 // Compute 2D f function
808 double f = 0.0;
809 double d0_f = 0.0;
810 double d1_f = 0.0;
811 for (int e = 0; e < 3; ++e) {
812 int i = e; // previous vertex
813 int j = (e + 1) % 3; // next vertex
814
815 // Get vertex information
816 std::array<double,3> &node_i = edgePoint_p[i];
817 std::array<double,3> &d0_node_i = d0_edgePoint_p[i];
818 std::array<double,3> &d1_node_i = d1_edgePoint_p[i];
819
820 std::array<double,3> &node_j = edgePoint_p[j];
821 std::array<double, 3> &d0_node_j = d0_edgeNormal_p[j];
822 std::array<double, 3> &d1_node_j = d1_edgeNormal_p[j];
823
824 std::array<double, 3> &normal_i = edgeNormal_p[i];
825 std::array<double, 3> &d0_normal_i = d0_edgeNormal_p[i];
826 std::array<double, 3> &d1_normal_i = d1_edgeNormal_p[i];
827
828 std::array<double, 3> &normal_j = edgeNormal_p[j];
829 std::array<double, 3> &d0_normal_j = d0_edgeNormal_p[j];
830 std::array<double, 3> &d1_normal_j = d1_edgeNormal_p[j];
831
832 // Compute lambda parameters
833 std::array<double, 3> edge = node_j - node_i;
834 std::array<double, 3> d0_edge = d0_node_j - d0_node_i;
835 std::array<double, 3> d1_edge = d1_node_j - d1_node_i;
836
837 double lambda_ij = -dotProduct(edge, normal_i) / dotProduct(normal_s, normal_i);
838 double d0_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d0_normal_i) + dotProduct(d0_edge, normal_i)) / dotProduct(normal_s, normal_i);
839 double d1_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d1_normal_i) + dotProduct(d1_edge, normal_i)) / dotProduct(normal_s, normal_i);
840
841 double lambda_ji = dotProduct(edge, normal_j) / dotProduct(normal_s, normal_j);
842 double d0_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d0_normal_j) + dotProduct(d0_edge, normal_j)) / dotProduct(normal_s, normal_j);
843 double d1_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d1_normal_j) + dotProduct(d1_edge, normal_j)) / dotProduct(normal_s, normal_j);
844
845 // Compute contribution to the f function
846 f += t[i] * t[j] * (lambda_ij * t[i] + lambda_ji * t[j]);
847 d0_f += (d0_t[i] * t[j] + t[i] * d0_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
848 + t[i] * t[j] * (d0_lambda_ij * t[i] + lambda_ij * d0_t[i] + d0_lambda_ji * t[j] + lambda_ji * d0_t[j]);
849 d1_f += (d1_t[i] * t[j] + t[i] * d1_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
850 + t[i] * t[j] * (d1_lambda_ij * t[i] + lambda_ij * d1_t[i] + d1_lambda_ji * t[j] + lambda_ji * d1_t[j]);
851 }
852
853 // Compute projection point on curved surface
854 (*projectionPoint) = point_m + f * normal_s;
855
856 // Compute projection normal on curved surface
857 std::array<double, 3> d0_point_p = d0_point_m + d0_f * normal_s;
858 std::array<double, 3> d1_point_p = d1_point_m + d1_f * normal_s;
859
860 (*projectionNormal) = crossProduct(d0_point_p, d1_point_p);
861 (*projectionNormal) /= norm2((*projectionNormal));
862
863 if (dotProduct(normal_s, (*projectionNormal)) < 0.0) {
864 (*projectionNormal) *= -1.0;
865 }
866}
867
889void LevelSetSegmentationSurfaceInfo::evalHighOrderProjectionOnPolygon(const std::array<double,3> &point,
890 const SegmentConstIterator &segmentItr,
891 std::array<double, 3> *projectionPoint,
892 std::array<double, 3> *projectionNormal) const
893{
894 // Get segment
895 const Cell &segment = *segmentItr;
896
897 // Get projection point on segment
898 ConstProxyVector<long> segmentVertexIds = m_surface->getFacetOrderedVertexIds(segment);
899
900 std::size_t nSegmentVertices = segmentVertexIds.size();
901 BITPIT_CREATE_WORKSPACE(segmentVertexCoords, std::array<double BITPIT_COMMA 3>, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
902 m_surface->getVertexCoords(segmentVertexIds.size(), segmentVertexIds.data(), segmentVertexCoords);
903
904 BITPIT_CREATE_WORKSPACE(tau, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
905 std::array<double, 3> point_s = CGElem::projectPointPolygon(point, nSegmentVertices, segmentVertexCoords, tau);
906
907 std::array<double, 3> lowOrderProjectionNormal = tau[0] * computeSegmentVertexNormal(segmentItr, 0, true);
908 for (std::size_t i = 1; i < nSegmentVertices; ++i) {
909 lowOrderProjectionNormal += tau[i] * computeSegmentVertexNormal(segmentItr, i, true);
910 }
911
912 // Early return if point_s consides with segment's node
913 double distanceTolerance = m_surface->getTol();
914 for (std::size_t i = 0; i < nSegmentVertices; ++i) {
915 if (utils::DoubleFloatingEqual()(tau[i], 1., distanceTolerance, distanceTolerance)) {
916 (*projectionPoint) = segmentVertexCoords[i];
917 (*projectionNormal) = computeSegmentVertexNormal(segmentItr, i, true);
918 return;
919 }
920 }
921
922 std::array<double, 3> d0_point_s;
923 std::array<double, 3> d1_point_s;
924
925 // Define parameterization of point_s
926 // Definition: point_s - x_central = t0 * (x_prev - x_central) + t1 * (x_next - x_central)
927 // where x_prev, x_central, x_next are subsequent polygon vertices
928 // Vertex x_central is the one corresponding to polygon angle closer to 90 degrees
929 double cosMin = std::numeric_limits<double>::max();
930 for (std::size_t i = 0; i < nSegmentVertices; ++i) {
931 int i_p = i; // previous
932 int i_c = (i + 1) % nSegmentVertices; // central
933 int i_n = (i + 2) % nSegmentVertices; // next
934
935 const std::array<double,3> &point_p = m_surface->getVertexCoords(segmentVertexIds[i_p]);
936 const std::array<double,3> &point_c = m_surface->getVertexCoords(segmentVertexIds[i_c]);
937 const std::array<double,3> &point_n = m_surface->getVertexCoords(segmentVertexIds[i_n]);
938
939 std::array<double,3> x_p = point_p - point_c;
940 std::array<double,3> x_n = point_n - point_c;
941
942 double cos = std::abs(dotProduct(x_p, x_n)) / (norm2(x_p) * norm2(x_n));
943 if (cos < cosMin) {
944 cosMin = cos;
945 d0_point_s = x_p;
946 d1_point_s = x_n;
947 }
948 }
949
950 // Get normal on segment
951 std::array<double, 3> facetNormal = computeSegmentNormal(segmentItr);
952 std::array<double, 3> normal_s = point - point_s;
953 double distance = norm2(normal_s);
954 if (utils::DoubleFloatingEqual()(distance, 0., distanceTolerance, distanceTolerance)) {
955 normal_s = facetNormal;
956 } else {
957 normal_s /= distance;
958 if (dotProduct(normal_s, facetNormal) < 0.0) {
959 normal_s *= -1.0;
960 }
961 }
962
963 // Compute projection point and normal at each edge "e"
964 int nSegmentEdges = nSegmentVertices;
965 BITPIT_CREATE_WORKSPACE(edgePoint_s, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
966 BITPIT_CREATE_WORKSPACE(d0_edgePoint_s, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
967 BITPIT_CREATE_WORKSPACE(d1_edgePoint_s, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
968 BITPIT_CREATE_WORKSPACE(edgePoint_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
969 BITPIT_CREATE_WORKSPACE(d0_edgePoint_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
970 BITPIT_CREATE_WORKSPACE(d1_edgePoint_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
971 BITPIT_CREATE_WORKSPACE(edgeNormal_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
972 BITPIT_CREATE_WORKSPACE(d0_edgeNormal_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
973 BITPIT_CREATE_WORKSPACE(d1_edgeNormal_p, std::array<double BITPIT_COMMA 3>, nSegmentEdges, ReferenceElementInfo::MAX_ELEM_EDGES);
974 for (int e = 0; e < nSegmentEdges; ++e) {
975 // Get edge's local and global node ids
976 ConstProxyVector<int> edgeLocalVertexIds = segment.getFaceLocalVertexIds(e);
977 int localId0 = edgeLocalVertexIds[0];
978 int localId1 = edgeLocalVertexIds[1];
979 int id0 = segment.getVertexId(localId0);
980 int id1 = segment.getVertexId(localId1);
981
982 // Get vertex information
983 const std::array<double,3> &node0 = m_surface->getVertexCoords(id0);
984 const std::array<double,3> &node1 = m_surface->getVertexCoords(id1);
985
986 std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, localId0, true);
987 std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, localId1, true);
988
989 // Get edge normal
990 std::array<double, 3> edgeNormal_s = computeSegmentEdgeNormal(segmentItr, e, true);
991
992 // Compute projection coefficients
993 std::array<double, 3> edge = node1 - node0;
994 double p01 = dotProduct(point_s - node0, edge) / dotProduct(edge, edge);
995 double d0_p01;
996 double d1_p01;
997 if (p01 > 1.0) {
998 p01 = 1.0;
999 d0_p01 = 0.0;
1000 d1_p01 = 0.0;
1001 } else if (p01 < 0.0) {
1002 p01 = 0.0;
1003 d0_p01 = 0.0;
1004 d1_p01 = 0.0;
1005 } else {
1006 d0_p01 = dotProduct(d0_point_s - node0, edge) / dotProduct(edge, edge);
1007 d1_p01 = dotProduct(d1_point_s - node0, edge) / dotProduct(edge, edge);
1008 }
1009 double p10 = 1.0 - p01;
1010 double d0_p10 = - d0_p01;
1011 double d1_p10 = - d1_p01;
1012
1013 // Compute lambda coefficients
1014 double lambda01 = -dotProduct(edge, normal0) / dotProduct(edgeNormal_s, normal0);
1015 double lambda10 = dotProduct(edge, normal1) / dotProduct(edgeNormal_s, normal1);
1016
1017 // Edge projection point
1018 double product1 = p10 * p01;
1019 double d_product1 = p01 - p10; // derivative wrt p10
1020
1021 double product2 = lambda01 * p10 + lambda10 * p01;
1022 double d_product2 = lambda01 - lambda10;
1023
1024 double f01 = product1 * product2;
1025 double d_f01 = d_product1 * product2 + product1 * d_product2;
1026 double d_d_f01 = 2.0 * (d_product1 * d_product2 - product2);
1027
1028 edgePoint_s[e] = p10 * node0 + p01 * node1;
1029 std::array<double, 3> d_edgePoint_s = node0 - node1;
1030
1031 d0_edgePoint_s[e] = d_edgePoint_s * d0_p10;
1032 d1_edgePoint_s[e] = d_edgePoint_s * d1_p10;
1033
1034 edgePoint_p[e] = edgePoint_s[e] + f01 * edgeNormal_s;
1035 std::array<double, 3> d_edgePoint_p = d_edgePoint_s + d_f01 * edgeNormal_s;
1036
1037 d0_edgePoint_p[e] = d_edgePoint_p * d0_p10;
1038 d1_edgePoint_p[e] = d_edgePoint_p * d1_p10;
1039
1040 // Edge surface normal
1041 double length = norm2(edge);
1042 edgeNormal_p[e] = d_f01 * edge + length * length * edgeNormal_s;
1043 std::array<double, 3> d_edgeNormal_p = d_d_f01 * edge;
1044
1045 double norm = norm2(edgeNormal_p[e]);
1046 double d_norm = dotProduct(edgeNormal_p[e], d_edgeNormal_p) / norm;
1047 edgeNormal_p[e] /= norm;
1048 d_edgeNormal_p = (d_edgeNormal_p - edgeNormal_p[e] * d_norm) / norm;
1049
1050 d0_edgeNormal_p[e] = d_edgeNormal_p * d0_p10;
1051 d1_edgeNormal_p[e] = d_edgeNormal_p * d1_p10;
1052 }
1053
1054 // Early return if point_s coinsides with a polygon's
1055 // edge (or a node of the intermediate polygon)
1056 for (int e = 0; e < nSegmentEdges; ++e) {
1057 double distance = norm2(edgePoint_s[e] - point_s);
1058 if (utils::DoubleFloatingEqual()(distance, 0.0, distanceTolerance, distanceTolerance)) {
1059 (*projectionPoint) = edgePoint_p[e];
1060 (*projectionNormal) = edgeNormal_p[e];
1061 return;
1062 }
1063 }
1064
1065 // Project surface point to the intermediator polygon
1066 // The projection is done by following formula:
1067 // point_m - point_s = sum i=0,vertices-1 { t[i] * (edgePoint_p[i] - edgePoint_s[i]) }
1068 // After various numerical experiments it is concluded that a quite smooth interpolation is
1069 // resulted from the following choice of the interpolation weights "t":
1070 // t[i] = w[i] / sum j=0,vertices-1 { w[j] }, where
1071 // w[i] = product j=0,vertices-1, j!=i { d[j]^4 }, where
1072 // d[i] is the distance between point_s and edgePoint_s[i]
1073 BITPIT_CREATE_WORKSPACE(d, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1074 BITPIT_CREATE_WORKSPACE(d0_d, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1075 BITPIT_CREATE_WORKSPACE(d1_d, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1076 for (std::size_t i = 0; i < nSegmentVertices; ++i) {
1077 d[i] = (edgePoint_s[i][0] - point_s[0]) * (edgePoint_s[i][0] - point_s[0])
1078 + (edgePoint_s[i][1] - point_s[1]) * (edgePoint_s[i][1] - point_s[1])
1079 + (edgePoint_s[i][2] - point_s[2]) * (edgePoint_s[i][2] - point_s[2]);
1080
1081 d[i] = d[i] * d[i];
1082
1083 d0_d[i] = 4.0 * d[i] * dotProduct((edgePoint_s[i] - point_s), (d0_edgePoint_s[i] - d0_point_s));
1084 d1_d[i] = 4.0 * d[i] * dotProduct((edgePoint_s[i] - point_s), (d1_edgePoint_s[i] - d1_point_s));
1085 }
1086
1087 BITPIT_CREATE_WORKSPACE(t, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1088 BITPIT_CREATE_WORKSPACE(d0_t, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1089 BITPIT_CREATE_WORKSPACE(d1_t, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1090 double sum = 0.0;
1091 double d0_sum = 0.0;
1092 double d1_sum = 0.0;
1093 for (std::size_t i = 0; i < nSegmentVertices; ++i) {
1094 t[i] = 1.0;
1095 d0_t[i] = 0.0;
1096 d1_t[i] = 0.0;
1097 for (std::size_t j = 0; j < nSegmentVertices; ++j) {
1098 if (i != j) {
1099 t[i] *= d[j];
1100
1101 double product = 1.0;
1102 for (std::size_t k = 0; k < nSegmentVertices; ++k) {
1103 if (j != k) {
1104 product *= d[k];
1105 }
1106 }
1107 d0_t[i] += d0_d[j] * product;
1108 d1_t[i] += d1_d[j] * product;
1109 }
1110 }
1111 sum += t[i];
1112 d0_sum += d0_t[i];
1113 d1_sum += d1_t[i];
1114 }
1115 for (std::size_t i = 0; i < nSegmentVertices; ++i) {
1116 t[i] /= sum;
1117 d0_t[i] = (d0_t[i] - t[i] * d0_sum) / sum;
1118 d1_t[i] = (d1_t[i] - t[i] * d1_sum) / sum;
1119 }
1120
1121 std::array<double, 3> point_m = {0.0, 0.0, 0.0};
1122 std::array<double, 3> d0_point_m = {0.0, 0.0, 0.0};
1123 std::array<double, 3> d1_point_m = {0.0, 0.0, 0.0};
1124 for (std::size_t i = 0; i < nSegmentVertices; ++i) {
1125 point_m += (edgePoint_p[i] - edgePoint_s[i]) * t[i];
1126 d0_point_m += d0_t[i] * (edgePoint_p[i]- edgePoint_s[i]) + t[i] * (d0_edgePoint_p[i] - d0_edgePoint_s[i]);
1127 d1_point_m += d1_t[i] * (edgePoint_p[i]- edgePoint_s[i]) + t[i] * (d1_edgePoint_p[i] - d1_edgePoint_s[i]);
1128 }
1129 point_m += point_s;
1130 d0_point_m += d0_point_s;
1131 d1_point_m += d1_point_s;
1132
1133 // Compute 2D f function
1134 double f = 0.0;
1135 double d0_f = 0.0;
1136 double d1_f = 0.0;
1137 for (int e = 0; e < nSegmentEdges; ++e) {
1138 int i = e; // previous vertex
1139 int j = (e + 1) % nSegmentEdges; // next vertex
1140
1141 // Get vertex information
1142 std::array<double,3> &node_i = edgePoint_p[i];
1143 std::array<double,3> &d0_node_i = d0_edgePoint_p[i];
1144 std::array<double,3> &d1_node_i = d1_edgePoint_p[i];
1145
1146 std::array<double, 3> &node_j = edgePoint_p[j];
1147 std::array<double, 3> &d0_node_j = d0_edgeNormal_p[j];
1148 std::array<double, 3> &d1_node_j = d1_edgeNormal_p[j];
1149
1150 std::array<double, 3> &normal_i = edgeNormal_p[i];
1151 std::array<double, 3> &d0_normal_i = d0_edgeNormal_p[i];
1152 std::array<double, 3> &d1_normal_i = d1_edgeNormal_p[i];
1153
1154 std::array<double, 3> &normal_j = edgeNormal_p[j];
1155 std::array<double, 3> &d0_normal_j = d0_edgeNormal_p[j];
1156 std::array<double, 3> &d1_normal_j = d1_edgeNormal_p[j];
1157
1158 // Compute lambda parameters
1159 std::array<double, 3> edge = node_j - node_i;
1160 std::array<double, 3> d0_edge = d0_node_j - d0_node_i;
1161 std::array<double, 3> d1_edge = d1_node_j - d1_node_i;
1162
1163 double lambda_ij = -dotProduct(edge, normal_i) / dotProduct(normal_s, normal_i);
1164 double d0_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d0_normal_i) + dotProduct(d0_edge, normal_i)) / dotProduct(normal_s, normal_i);
1165 double d1_lambda_ij = -(dotProduct(edge - lambda_ij * normal_s, d1_normal_i) + dotProduct(d1_edge, normal_i)) / dotProduct(normal_s, normal_i);
1166
1167 double lambda_ji = dotProduct(edge, normal_j) / dotProduct(normal_s, normal_j);
1168 double d0_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d0_normal_j) + dotProduct(d0_edge, normal_j)) / dotProduct(normal_s, normal_j);
1169 double d1_lambda_ji = (dotProduct(edge - lambda_ji * normal_s, d1_normal_j) + dotProduct(d1_edge, normal_j)) / dotProduct(normal_s, normal_j);
1170
1171 // Compute contribution to the f function
1172 f += t[i] * t[j] * (lambda_ij * t[i] + lambda_ji * t[j]);
1173 d0_f += (d0_t[i] * t[j] + t[i] * d0_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
1174 + t[i] * t[j] * (d0_lambda_ij * t[i] + lambda_ij * d0_t[i] + d0_lambda_ji * t[j] + lambda_ji * d0_t[j]);
1175 d1_f += (d1_t[i] * t[j] + t[i] * d1_t[j]) * (lambda_ij * t[i] + lambda_ji * t[j])
1176 + t[i] * t[j] * (d1_lambda_ij * t[i] + lambda_ij * d1_t[i] + d1_lambda_ji * t[j] + lambda_ji * d1_t[j]);
1177 }
1178
1179 // Compute projection point on curved surface
1180 (*projectionPoint) = point_m + f * normal_s;
1181
1182 // Compute projection normal on curved surface
1183 std::array<double, 3> d0_point_p = d0_point_m + d0_f * normal_s;
1184 std::array<double, 3> d1_point_p = d1_point_m + d1_f * normal_s;
1185
1186 (*projectionNormal) = crossProduct(d0_point_p, d1_point_p);
1187 (*projectionNormal) /= norm2((*projectionNormal));
1188
1189 if (dotProduct(normal_s, (*projectionNormal)) < 0.0) {
1190 (*projectionNormal) *= -1.0;
1191 }
1192}
1193
1213void LevelSetSegmentationSurfaceInfo::evalHighOrderProjection(const std::array<double,3> &point,
1214 const SegmentConstIterator &segmentItr,
1215 std::array<double, 3> *projectionPoint,
1216 std::array<double, 3> *projectionNormal) const
1217{
1218
1219 const Cell &segment = *segmentItr;
1220 ElementType segmentType = segment.getType();
1221 switch (segmentType) {
1222
1223 case ElementType::VERTEX:
1224 {
1225 evalProjectionOnVertex(point, segmentItr, projectionPoint, projectionNormal);
1226 return;
1227 }
1228
1229 case ElementType::LINE:
1230 {
1231 evalHighOrderProjectionOnLine(point, segmentItr, projectionPoint, projectionNormal);
1232 return;
1233 }
1234
1235 case ElementType::TRIANGLE:
1236 {
1237 evalHighOrderProjectionOnTriangle(point, segmentItr, projectionPoint, projectionNormal);
1238 return;
1239 }
1240
1241 default:
1242 {
1243 evalHighOrderProjectionOnPolygon(point, segmentItr, projectionPoint, projectionNormal);
1244 return;
1245 }
1246
1247 }
1248}
1249
1258void LevelSetSegmentationSurfaceInfo::evalLowOrderProjectionOnLine(const std::array<double, 3> &point,
1259 const SegmentConstIterator &segmentItr,
1260 std::array<double, 3> *projectionPoint,
1261 std::array<double, 3> *projectionNormal) const
1262{
1263 const Cell &segment = *segmentItr;
1264 assert(segment.getType() == ElementType::LINE);
1265
1266 ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
1267
1268 const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
1269 const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
1270
1271 int nSegmentVertices = segment.getVertexCount();
1272 BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1273 (*projectionPoint) = CGElem::projectPointSegment(point, point0, point1, lambda);
1274
1275 std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, 0, true);
1276 std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, 1, true);
1277
1278 (*projectionNormal) = lambda[0] * normal0 + lambda[1] * normal1;
1279 (*projectionNormal) /= norm2((*projectionNormal));
1280}
1281
1290void LevelSetSegmentationSurfaceInfo::evalLowOrderProjectionOnTriangle(const std::array<double, 3> &point,
1291 const SegmentConstIterator &segmentItr,
1292 std::array<double, 3> *projectionPoint,
1293 std::array<double, 3> *projectionNormal) const
1294{
1295 const Cell &segment = *segmentItr;
1296 assert(segment.getType() == ElementType::TRIANGLE);
1297
1298 ConstProxyVector<long> segmentVertexIds = segment.getVertexIds();
1299
1300 const std::array<double,3> &point0 = m_surface->getVertexCoords(segmentVertexIds[0]);
1301 const std::array<double,3> &point1 = m_surface->getVertexCoords(segmentVertexIds[1]);
1302 const std::array<double,3> &point2 = m_surface->getVertexCoords(segmentVertexIds[2]);
1303
1304 int nSegmentVertices = segment.getVertexCount();
1305 BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1306 (*projectionPoint) = CGElem::projectPointTriangle(point, point0, point1, point2, lambda);
1307
1308 std::array<double, 3> normal0 = computeSegmentVertexNormal(segmentItr, 0, true);
1309 std::array<double, 3> normal1 = computeSegmentVertexNormal(segmentItr, 1, true);
1310 std::array<double, 3> normal2 = computeSegmentVertexNormal(segmentItr, 2, true);
1311
1312 (*projectionNormal) = lambda[0] * normal0 + lambda[1] * normal1 + lambda[2] * normal2;
1313 (*projectionNormal) /= norm2((*projectionNormal));
1314}
1315
1324void LevelSetSegmentationSurfaceInfo::evalLowOrderProjectionOnPolygon(const std::array<double, 3> &point,
1325 const SegmentConstIterator &segmentItr,
1326 std::array<double, 3> *projectionPoint,
1327 std::array<double, 3> *projectionNormal) const
1328{
1329 const Cell &segment = *segmentItr;
1330
1331 ConstProxyVector<long> segmentVertexIds = m_surface->getFacetOrderedVertexIds(segment);
1332
1333 std::size_t nSegmentVertices = segmentVertexIds.size();
1334 BITPIT_CREATE_WORKSPACE(segmentVertexCoords, std::array<double BITPIT_COMMA 3>, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1335 m_surface->getVertexCoords(segmentVertexIds.size(), segmentVertexIds.data(), segmentVertexCoords);
1336
1337 BITPIT_CREATE_WORKSPACE(lambda, double, nSegmentVertices, ReferenceElementInfo::MAX_ELEM_VERTICES);
1338 (*projectionPoint) = CGElem::projectPointPolygon(point, nSegmentVertices, segmentVertexCoords, lambda);
1339
1340 (*projectionNormal) = lambda[0] * computeSegmentVertexNormal(segmentItr, 0, true);
1341 for (std::size_t i = 1; i < nSegmentVertices; ++i) {
1342 (*projectionNormal) += lambda[i] * computeSegmentVertexNormal(segmentItr, i, true);
1343 }
1344 (*projectionNormal) /= norm2(*projectionNormal);
1345}
1346
1355void LevelSetSegmentationSurfaceInfo::evalLowOrderProjection(const std::array<double, 3> &point,
1356 const SegmentConstIterator &segmentItr,
1357 std::array<double, 3> *projectionPoint,
1358 std::array<double, 3> *projectionNormal) const
1359{
1360 const Cell &segment = *segmentItr;
1361 ElementType segmentType = segment.getType();
1362 switch (segmentType) {
1363
1364 case ElementType::VERTEX:
1365 {
1366 evalProjectionOnVertex(point, segmentItr, projectionPoint, projectionNormal);
1367 return;
1368 }
1369
1370 case ElementType::LINE:
1371 {
1372 evalLowOrderProjectionOnLine(point, segmentItr, projectionPoint, projectionNormal);
1373 return;
1374 }
1375
1376 case ElementType::TRIANGLE:
1377 {
1378 evalLowOrderProjectionOnTriangle(point, segmentItr, projectionPoint, projectionNormal);
1379 return;
1380 }
1381
1382 default:
1383 {
1384 evalLowOrderProjectionOnPolygon(point, segmentItr, projectionPoint, projectionNormal);
1385 return;
1386 }
1387
1388 }
1389}
1390
1399void LevelSetSegmentationSurfaceInfo::evalProjection(const std::array<double, 3> &point,
1400 const SegmentConstIterator &segmentItr,
1401 std::array<double, 3> *projectionPoint,
1402 std::array<double, 3> *projectionNormal) const
1403{
1404 if (m_surfaceSmoothing == LevelSetSurfaceSmoothing::HIGH_ORDER) {
1405 evalHighOrderProjection(point, segmentItr, projectionPoint, projectionNormal);
1406 } else {
1407 evalLowOrderProjection(point, segmentItr, projectionPoint, projectionNormal);
1408 }
1409}
1410
1418void LevelSetSegmentationSurfaceInfo::evalProjection(const std::array<double, 3> &point,
1419 const SegmentConstIterator &segmentItr,
1420 std::array<double, 3> *projectionPoint) const
1421{
1422 std::array<double, 3> projectionNormal;
1423 evalProjection(point, segmentItr, projectionPoint, &projectionNormal);
1424
1425 BITPIT_UNUSED(projectionNormal);
1426}
1427
1436std::array<double, 3> LevelSetSegmentationSurfaceInfo::evalProjection(const std::array<double, 3> &point,
1437 const SegmentConstIterator &segmentItr,
1438 double *lambda) const
1439{
1440 std::array<double, 3> projectionPoint;
1441 evalProjection(point, segmentItr, &projectionPoint);
1442
1443 const Cell &segment = *segmentItr;
1444 m_surface->evalBarycentricCoordinates(segment.getId(), point, lambda);
1445 return projectionPoint;
1446}
1447
1470std::array<double,3> LevelSetSegmentationSurfaceInfo::computePseudoNormal(const SegmentConstIterator &segmentItr,
1471 const double *lambda ) const
1472{
1473 // Early return if the segment is a point
1474 const Cell &segment = *segmentItr;
1475 ElementType segmentType = segment.getType();
1476 if (segmentType == ElementType::VERTEX) {
1477 return {{0., 0., 0.}};
1478 }
1479
1480 // Evaluate pseudo normal
1481 int positionFlag;
1482 if (segmentType == ElementType::LINE) {
1483 positionFlag = CGElem::convertBarycentricToFlagSegment(lambda, m_surface->getTol());
1484 } else {
1485 int nSegmentVertices = segment.getVertexCount();
1486 positionFlag = CGElem::convertBarycentricToFlagPolygon(nSegmentVertices, lambda, m_surface->getTol());
1487 if (positionFlag > 0) {
1488 int polygonVertex = positionFlag - 1;
1489 int elementVertex = m_surface->getFacetOrderedLocalVertex(segment, polygonVertex);
1490
1491 positionFlag = elementVertex + 1;
1492 } else if (positionFlag < 0) {
1493 int polygonEdge = (- positionFlag) - 1;
1494 int elementEdge = m_surface->getFacetOrderedLocalEdge(segment, polygonEdge);
1495
1496 positionFlag = - (elementEdge + 1);
1497 }
1498 }
1499
1500 std::array<double,3> pseudoNormal;
1501 if (positionFlag == 0) {
1502 pseudoNormal = computeSegmentNormal(segmentItr);
1503 } else if (positionFlag > 0) {
1504 int vertex = positionFlag - 1;
1505 pseudoNormal = computeSegmentVertexNormal(segmentItr, vertex, false);
1506 } else {
1507 int edge = (- positionFlag) - 1;
1508 pseudoNormal = computeSegmentEdgeNormal(segmentItr, edge, false);
1509 }
1510
1511 return pseudoNormal;
1512}
1513
1526std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSurfaceNormal(const SegmentConstIterator &segmentItr,
1527 const double *lambda ) const
1528{
1529 // Early return if the segment is a point
1530 const Cell &segment = *segmentItr;
1531 ElementType segmentType = segment.getType();
1532 if (segmentType == ElementType::VERTEX) {
1533 return {{0., 0., 0.}};
1534 }
1535
1536 // Evaluate surface normal
1537 std::size_t nSegmentVertices = segment.getVertexCount();
1538 std::array<double,3> surfaceNormal = lambda[0] * computeSegmentVertexNormal(segmentItr, 0, true);
1539 for (std::size_t i = 1; i < nSegmentVertices; ++i) {
1540 surfaceNormal += lambda[i] * computeSegmentVertexNormal(segmentItr, i, true);
1541 }
1542 surfaceNormal /= norm2(surfaceNormal);
1543
1544 return surfaceNormal;
1545}
1546
1555std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSegmentNormal(const SegmentConstIterator &segmentItr ) const {
1556
1557 std::size_t segmentRawId = segmentItr.getRawIndex();
1558 std::array<double, 3> *segmentNormal = m_segmentNormalsStorage.rawData(segmentRawId);
1559 if (!m_segmentNormalsValid.rawAt(segmentRawId)) {
1560 *segmentNormal = m_surface->evalFacetNormal(segmentItr->getId());
1561 m_segmentNormalsValid.rawAt(segmentRawId) = true;
1562 }
1563
1564 return *segmentNormal;
1565}
1566
1576std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSegmentEdgeNormal(const SegmentConstIterator &segmentItr, int edge, bool limited ) const {
1577
1578 long segmentId = segmentItr.getId();
1579
1580 std::array<double, 3> limitedEdgeNormal;
1581 std::array<double, 3> unlimitedEdgeNormal;
1582 m_surface->evalEdgeNormals(segmentId, edge, m_featureAngle, &unlimitedEdgeNormal, &limitedEdgeNormal) ;
1583
1584 if (limited) {
1585 return limitedEdgeNormal;
1586 }
1587
1588 return unlimitedEdgeNormal;
1589}
1590
1602std::array<double,3> LevelSetSegmentationSurfaceInfo::computeSegmentVertexNormal(const SegmentConstIterator &segmentItr, int vertex, bool limited ) const {
1603
1604 // Segment information
1605 long segmentId = segmentItr.getId();
1606 long segmentRawId = segmentItr.getRawIndex();
1607 const Cell &segment = *segmentItr;
1608
1609 // Update the cache
1610 //
1611 // Cache is updated for all the vertices of the segment.
1612 long vertexId = segment.getVertexId(vertex);
1613 std::size_t vertexRawId = m_surface->getVertices().getRawIndex(vertexId);
1614 bool hasUnlimitedNormal = m_unlimitedVertexNormalsValid.rawAt(vertexRawId);
1615
1616 bool hasLimitedNormal = m_limitedSegmentVertexNormalValid[m_segmentVertexOffset.rawAt(segmentRawId) + vertex];
1617
1618 if (!hasUnlimitedNormal || !hasLimitedNormal) {
1619 static std::vector<long> vertexNeighbours;
1620 vertexNeighbours.clear();
1621 m_surface->findCellVertexNeighs(segmentId, vertex, &vertexNeighbours);
1622
1623 std::array<double, 3> limitedVertexNormal;
1624 std::array<double, 3> unlimitedVertexNormal;
1625 if (hasUnlimitedNormal) {
1626 limitedVertexNormal = m_surface->evalLimitedVertexNormal(segmentId, vertex, vertexNeighbours.size(), vertexNeighbours.data(), m_featureAngle) ;
1627 unlimitedVertexNormal = m_unlimitedVertexNormalsStorage.rawAt(vertexRawId);
1628 } else {
1629 m_surface->evalVertexNormals(segmentId, vertex, vertexNeighbours.size(), vertexNeighbours.data(), m_featureAngle, &unlimitedVertexNormal, &limitedVertexNormal) ;
1630 }
1631
1632 // Store vertex limited normal
1633 //
1634 // Both limited and unlimited normal are evaluated, however limited
1635 // normal is only stored if its misalignment with respect to the
1636 // unlimited normal is greater than a defined tolerance.
1637 if (!hasLimitedNormal) {
1638 double misalignment = norm2(unlimitedVertexNormal - limitedVertexNormal) ;
1639 if (misalignment >= m_surface->getTol()) {
1640 std::pair<long, int> segmentVertexKey = std::make_pair(segmentId, vertex);
1641 m_limitedSegmentVertexNormalStorage.insert({segmentVertexKey, std::move(limitedVertexNormal)}) ;
1642 }
1643 m_limitedSegmentVertexNormalValid[m_segmentVertexOffset.rawAt(segmentRawId) + vertex] = true;
1644 }
1645
1646 // Store vertex unlimited normal
1647 if (!hasUnlimitedNormal ) {
1648 m_unlimitedVertexNormalsStorage.rawAt(vertexRawId) = std::move(unlimitedVertexNormal);
1649 m_unlimitedVertexNormalsValid.rawAt(vertexRawId) = true ;
1650 }
1651 }
1652
1653 // Get the normal from the cache
1654 if (limited) {
1655 std::pair<long, int> segmentVertexKey = std::make_pair(segmentId, vertex);
1656 auto limitedSegmentVertexNormal = m_limitedSegmentVertexNormalStorage.find(segmentVertexKey);
1657 if (limitedSegmentVertexNormal != m_limitedSegmentVertexNormalStorage.end()) {
1658 return limitedSegmentVertexNormal->second;
1659 }
1660 }
1661
1662 return m_unlimitedVertexNormalsStorage.rawAt(vertexRawId);
1663}
1664
1677const double LevelSetSegmentationBaseObject::AUTOMATIC_SEARCH_RADIUS = -1;
1678
1687std::size_t LevelSetSegmentationBaseObject::createFieldCellCache(LevelSetField field, std::size_t cacheId)
1688{
1689 switch(field) {
1690
1691 case LevelSetField::SUPPORT:
1692 return createFieldCellCache<long>(field, cacheId);
1693
1694 default:
1695 return LevelSetObject::createFieldCellCache(field, cacheId);
1696
1697 }
1698}
1699
1704LevelSetFieldset LevelSetSegmentationBaseObject::getSupportedFields() const
1705{
1706 LevelSetFieldset supportedFields = LevelSetObject::getSupportedFields();
1707 supportedFields.push_back(LevelSetField::PART);
1708 supportedFields.push_back(LevelSetField::NORMAL);
1709 supportedFields.push_back(LevelSetField::SUPPORT);
1710
1711 return supportedFields;
1712}
1713
1720void LevelSetSegmentationBaseObject::fillFieldCellCache(LevelSetField field, long id)
1721{
1722 switch (field) {
1723
1724 case LevelSetField::SUPPORT:
1725 evalCellSupport(id);
1726 break;
1727
1728 default:
1729 LevelSetObject::fillFieldCellCache(field, id);
1730
1731 }
1732}
1733
1740const SurfUnstructured & LevelSetSegmentationBaseObject::evalCellSurface(long id) const
1741{
1742 return _evalCellSurface(id);
1743}
1744
1751int LevelSetSegmentationBaseObject::evalCellPart(long id) const
1752{
1753 auto evaluator = [this] (long id)
1754 {
1755 return _evalCellPart(id);
1756 };
1757
1758 auto fallback = [] (long id)
1759 {
1760 BITPIT_UNUSED(id);
1761
1762 return levelSetDefaults::PART;
1763 };
1764
1765 LevelSetField field = LevelSetField::PART;
1766 int part = evalCellFieldCached<int>(field, id, evaluator, fallback);
1767
1768 return part;
1769}
1770
1806LevelSetIntersectionStatus LevelSetSegmentationBaseObject::_intersectSurface(long id, double distance, LevelSetIntersectionMode mode) const
1807{
1808 // Get surface information
1809 const SurfUnstructured &surface = evalCellSurface(id);
1810
1811 // Early return if the surface is empty
1812 if (surface.empty()) {
1813 return LevelSetIntersectionStatus::FALSE;
1814 }
1815
1816 // Evaluate intersection using base class
1817 return LevelSetObject::_intersectSurface(id, distance, mode);
1818}
1819
1827std::array<double,3> LevelSetSegmentationBaseObject::evalCellNormal(long id, bool signedLevelSet) const
1828{
1829 // Evaluate signed normal
1830 //
1831 // The normal stored in the cache is unsigned.
1832 auto evaluator = [this] (long id)
1833 {
1834 return _evalCellNormal(id, false);
1835 };
1836
1837 auto fallback = [] (long id)
1838 {
1839 BITPIT_UNUSED(id);
1840
1841 return levelSetDefaults::NORMAL;
1842 };
1843
1844 LevelSetField field = LevelSetField::NORMAL;
1845 std::array<double, 3> normal = evalCellFieldCached<std::array<double, 3>>(field, id, evaluator, fallback);
1846
1847 // Evaluate the normal with the correct signdness
1848 //
1849 // If an unsigned evaluation is requested, the orientation of the surface should be discarded
1850 // and in order to have a normal that is agnostic with respect the two sides of the surface.
1851 if (signedLevelSet) {
1852 short cellSign = evalCellSign(id);
1853 if (cellSign <= 0) {
1854 normal *= static_cast<double>(cellSign);;
1855 }
1856 }
1857
1858 return normal;
1859}
1860
1875long LevelSetSegmentationBaseObject::evalCellSupport(long id, double searchRadius) const
1876{
1877 // Evaluate support
1878 //
1879 // If we are using a limited search radius, it is not possible to use the support stored
1880 // in the cache, because it may be outside the search radius.
1881 auto evaluator = [this, searchRadius] (long id)
1882 {
1883 // Evaluate the search radius for support evaluation
1884 //
1885 // For the automatic evaluation of the search range, it is possible to use the
1886 // information about the cell zone:
1887 // - cells that are inside the narrow band because their distance from the surface is
1888 // less than the narrow band size can use a search radius equal to the narrow band
1889 // size
1890 // - cells that are inside the narrow band because they intersects the surface can use
1891 // a search radius equal to their bounding radius;
1892 // - cells that are inside the narrow band because their one of their face neighbours
1893 // intersect the surface can use a search radius equal to their bounding radius plus
1894 // the bounding diameter of the smallest neighbour that intersects the surface.
1895 double supportSearchRadius;
1896 if (searchRadius == AUTOMATIC_SEARCH_RADIUS) {
1897 LevelSetCellLocation cellLocation = getCellLocation(id);
1898 if (cellLocation == LevelSetCellLocation::NARROW_BAND_DISTANCE) {
1899 supportSearchRadius = this->m_narrowBandSize;
1900 } else if (cellLocation == LevelSetCellLocation::NARROW_BAND_INTERSECTED) {
1901 supportSearchRadius = m_kernel->computeCellBoundingRadius(id);
1902 } else if (cellLocation == LevelSetCellLocation::NARROW_BAND_NEIGHBOUR) {
1903 // Evaluate the bounding diameter of the smallest intersected face neighbour
1904 supportSearchRadius = std::numeric_limits<double>::max();
1905 auto neighProcessor = [this, &supportSearchRadius](long neighId, int layer) {
1906 BITPIT_UNUSED(layer);
1907
1908 // Discard neighbours that are not intersected
1909 LevelSetCellLocation neighLocation = getCellLocation(neighId);
1910 if (neighLocation != LevelSetCellLocation::NARROW_BAND_INTERSECTED) {
1911 return false;
1912 }
1913
1914 // Consider the bounding diameter of the smallest intersected neighbour
1915 supportSearchRadius = std::min(2 * m_kernel->computeCellBoundingRadius(neighId), supportSearchRadius);
1916
1917 // Continue processing the other neighbours
1918 return false;
1919 };
1920
1921 const VolumeKernel &mesh = *(m_kernel->getMesh()) ;
1922 mesh.processCellFaceNeighbours(id, 1, neighProcessor);
1923
1924 // Ad the bounding radius of the cell
1925 supportSearchRadius += m_kernel->computeCellBoundingRadius(id);
1926 } else {
1927 supportSearchRadius = std::numeric_limits<double>::max();
1928 }
1929 } else {
1930 supportSearchRadius = searchRadius;
1931 }
1932
1933 // Evaluate cell support
1934 return _evalCellSupport(id, supportSearchRadius);
1935 };
1936
1937 auto fallback = [] (long id)
1938 {
1939 BITPIT_UNUSED(id);
1940
1941 return levelSetDefaults::SUPPORT;
1942 };
1943
1944 LevelSetField field = LevelSetField::SUPPORT;
1945
1946 long support;
1947 if (searchRadius < std::numeric_limits<double>::max() && searchRadius != AUTOMATIC_SEARCH_RADIUS) {
1948 // Evaluate the support from scratch
1949 support = evalCellField<long>(field, id, evaluator, fallback);
1950
1951 // Update the cache if the support is valid
1952 if (support >= 0) {
1953 fillFieldCellCache(field, id, support);
1954 }
1955 } else {
1956 // Evaluate the support
1957 support = evalCellFieldCached<long>(field, id, evaluator, fallback);
1958 }
1959
1960 return support;
1961}
1962
1969const SurfUnstructured & LevelSetSegmentationBaseObject::evalSurface(const std::array<double,3> &point) const
1970{
1971 return _evalSurface(point);
1972}
1973
1980int LevelSetSegmentationBaseObject::evalPart(const std::array<double,3> &point) const
1981{
1982 return _evalPart(point);
1983}
1984
1992std::array<double,3> LevelSetSegmentationBaseObject::evalNormal(const std::array<double,3> &point, bool signedLevelSet) const
1993{
1994 // Project the point on the surface and evaluate the point-projection vector
1995 std::array<double, 3> projectionPoint;
1996 std::array<double, 3> projectionNormal;
1997 evalProjection(point, signedLevelSet, &projectionPoint, &projectionNormal);
1998
1999 return projectionNormal;
2000}
2001
2008long LevelSetSegmentationBaseObject::evalSupport(const std::array<double,3> &point) const
2009{
2010 return _evalSupport(point);
2011}
2012
2021long LevelSetSegmentationBaseObject::evalSupport(const std::array<double,3> &point, double searchRadius) const
2022{
2023 return _evalSupport(point, searchRadius);
2024}
2025
2037void LevelSetSegmentationBaseObject::evalProjection(const std::array<double,3> &point,
2038 bool signedLevelSet,
2039 std::array<double, 3> *projectionPoint,
2040 std::array<double, 3> *projectionNormal) const
2041{
2042 _evalProjection(point, signedLevelSet, projectionPoint, projectionNormal);
2043}
2044
2051int LevelSetSegmentationBaseObject::_evalCellPart(long id) const
2052{
2053 long support = evalCellSupport(id);
2054 const SurfUnstructured &surface = evalCellSurface(id);
2055 int part = surface.getCell(support).getPID();
2056
2057 return part;
2058}
2059
2066int LevelSetSegmentationBaseObject::_evalPart(const std::array<double,3> &point) const
2067{
2068 long support = evalSupport(point);
2069 const SurfUnstructured &surface = evalSurface(point);
2070 int part = surface.getCell(support).getPID();
2071
2072 return part;
2073}
2074
2081void LevelSetSegmentationBaseObject::addVTKOutputData( LevelSetField field, const std::string &objectName)
2082{
2083 VTK &vtkWriter = this->m_kernel->getMesh()->getVTK() ;
2084 std::string name = this->getVTKOutputDataName(field, objectName);
2085
2086 switch (field) {
2087
2088 case LevelSetField::SUPPORT:
2089 vtkWriter.addData<long>( name, VTKFieldType::SCALAR, VTKLocation::CELL, this);
2090 break;
2091
2092 case LevelSetField::PART:
2093 vtkWriter.addData<int>( name, VTKFieldType::SCALAR, VTKLocation::CELL, this);
2094 break;
2095
2096 case LevelSetField::NORMAL:
2097 vtkWriter.addData<double>( name, VTKFieldType::VECTOR, VTKLocation::CELL, this);
2098 break;
2099
2100 default:
2101 LevelSetObject::addVTKOutputData(field, objectName);
2102 break;
2103
2104 }
2105}
2106
2113std::string LevelSetSegmentationBaseObject::getVTKOutputFieldName( LevelSetField field) const
2114{
2115 switch (field) {
2116
2117 case LevelSetField::SUPPORT:
2118 return "SupportId";
2119
2120 case LevelSetField::PART:
2121 return "PartId";
2122
2123 case LevelSetField::NORMAL:
2124 return "Normal";
2125
2126 default:
2127 return LevelSetObject::getVTKOutputFieldName(field);
2128
2129 }
2130}
2131
2141void LevelSetSegmentationBaseObject::flushVTKOutputData(std::fstream &stream, VTKFormat format,
2142 LevelSetField field) const
2143{
2144 switch (field) {
2145
2146 case LevelSetField::SUPPORT:
2147 {
2148 auto evaluator = [this] (long id) { return evalCellSupport(id); };
2149 auto fallback = [] (long id) { BITPIT_UNUSED(id); return levelSetDefaults::SUPPORT; };
2150 flushVTKOutputData<double>(stream, format, field, evaluator, fallback);
2151 break;
2152 }
2153
2154 case LevelSetField::PART:
2155 {
2156 auto evaluator = [this] (long id) { return evalCellPart(id); };
2157 auto fallback = [] (long id) { BITPIT_UNUSED(id); return levelSetDefaults::PART; };
2158 flushVTKOutputData<double>(stream, format, field, evaluator, fallback);
2159 break;
2160 }
2161
2162 case LevelSetField::NORMAL:
2163 {
2164 auto evaluator = [this] (long id) { return evalCellNormal(id, true); };
2165 auto fallback = [] (long id) { BITPIT_UNUSED(id); return levelSetDefaults::NORMAL; };
2166 flushVTKOutputData<double>(stream, format, field, evaluator, fallback);
2167 break;
2168 }
2169
2170 default:
2171 LevelSetObject::flushVTKOutputData(stream, format, field);
2172 break;
2173
2174 }
2175}
2176
2183int LevelSetSegmentationBaseObject::getPart(long cellId) const
2184{
2185 return evalCellPart(cellId);
2186}
2187
2194std::array<double,3> LevelSetSegmentationBaseObject::getNormal(long cellId) const
2195{
2196 return evalCellNormal(cellId, m_defaultSignedLevelSet);
2197}
2198
2205long LevelSetSegmentationBaseObject::getSupport(long cellId) const
2206{
2207 return evalCellSupport(cellId);
2208}
2209
2216double LevelSetSegmentationBaseObject::getSurfaceFeatureSize(long cellId) const
2217{
2218 long support = evalCellSupport(cellId);
2219 if (support < 0) {
2220 return (- levelSetDefaults::SIZE);
2221 }
2222
2223 const SurfUnstructured &surface = evalCellSurface(cellId);
2224
2225 return surface.evalCellSize(support);
2226}
2227
2234int LevelSetSegmentationBaseObject::getPart(const std::array<double,3> &point) const
2235{
2236 return evalPart(point);
2237}
2238
2245std::array<double,3> LevelSetSegmentationBaseObject::getNormal(const std::array<double,3> &point) const
2246{
2247 return evalNormal(point, m_defaultSignedLevelSet);
2248}
2249
2256long LevelSetSegmentationBaseObject::getSupport(const std::array<double,3> &point) const
2257{
2258 return evalSupport(point);
2259}
2260
2267double LevelSetSegmentationBaseObject::getSurfaceFeatureSize(const std::array<double,3> &point) const
2268{
2269 long support = evalSupport(point);
2270 if (support < 0) {
2271 return (- levelSetDefaults::SIZE);
2272 }
2273
2274 const SurfUnstructured &surface = evalSurface(point);
2275
2276 return surface.evalCellSize(support);
2277}
2278
2289LevelSetSegmentationObject::LevelSetSegmentationObject(int id)
2290 : LevelSetSegmentationBaseObject(id),
2291 m_surfaceInfo(nullptr)
2292{
2293}
2294
2301LevelSetSegmentationObject::LevelSetSegmentationObject(int id, std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle)
2302 : LevelSetSegmentationBaseObject(id)
2303{
2304 setSurface(std::move(surface), featureAngle);
2305}
2306
2313LevelSetSegmentationObject::LevelSetSegmentationObject(int id, const SurfUnstructured *surface, double featureAngle)
2314 : LevelSetSegmentationBaseObject(id)
2315{
2316 setSurface(surface, featureAngle);
2317}
2318
2326LevelSetSegmentationObject::LevelSetSegmentationObject(int id, const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing)
2327 : LevelSetSegmentationBaseObject(id)
2328{
2329 setSurface(surface, featureAngle, surfaceSmoothing);
2330}
2331
2337LevelSetSegmentationObject::LevelSetSegmentationObject(const LevelSetSegmentationObject &other)
2338 : LevelSetSegmentationBaseObject(other)
2339{
2340 if (other.m_surfaceInfo) {
2341 m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(*(other.m_surfaceInfo)));
2342 } else {
2343 m_surfaceInfo = nullptr;
2344 }
2345}
2346
2352bool LevelSetSegmentationObject::empty() const
2353{
2354 return getSurface().empty();
2355}
2356
2364
2369const SurfUnstructured & LevelSetSegmentationObject::getSurface() const {
2370 return m_surfaceInfo->getSurface();
2371}
2372
2387void LevelSetSegmentationObject::setSurface(std::unique_ptr<const SurfUnstructured> &&surface, bool force){
2388 setSurface(std::move(surface), LevelSetSegmentationSurfaceInfo::DEFAULT_FEATURE_ANGLE, force);
2389}
2390
2404void LevelSetSegmentationObject::setSurface(std::unique_ptr<const SurfUnstructured> &&surface, double featureAngle, bool force){
2405 if (m_surfaceInfo) {
2406 // Check if replacing an existing surface is allowed
2407 if (!force) {
2408 throw std::runtime_error ("The surface can only be set once.");
2409 }
2410
2411 // Replace the surface
2412 m_surfaceInfo->setSurface(std::move(surface), featureAngle);
2413 } else {
2414 // Set surface
2415 //
2416 // Since this is the first time we set the surface, there is no need
2417 // to clear the caches.
2418 m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(std::move(surface), featureAngle));
2419 }
2420}
2421
2439
2453void LevelSetSegmentationObject::setSurface(const SurfUnstructured *surface, double featureAngle, bool force){
2454 if (m_surfaceInfo) {
2455 // Check if replacing an existing surface is allowed
2456 if (!force) {
2457 throw std::runtime_error ("The surface can only be set once.");
2458 }
2459
2460 // Replace the surface
2461 m_surfaceInfo->setSurface(surface, featureAngle);
2462 } else {
2463 // Set surface
2464 //
2465 // Since this is the first time we set the surface, there is no need
2466 // to clear the caches.
2467 m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(surface, featureAngle));
2468 }
2469}
2470
2485void LevelSetSegmentationObject::setSurface(const SurfUnstructured *surface, double featureAngle, LevelSetSurfaceSmoothing surfaceSmoothing, bool force){
2486 if (m_surfaceInfo) {
2487 // Check if replacing an existing surface is allowed
2488 if (!force) {
2489 throw std::runtime_error ("The surface can only be set once.");
2490 }
2491
2492 // Replace the surface
2493 m_surfaceInfo->setSurface(surface, featureAngle, surfaceSmoothing);
2494 } else {
2495 // Set surface
2496 //
2497 // Since this is the first time we set the surface, there is no need
2498 // to clear the caches.
2499 m_surfaceInfo = std::unique_ptr<LevelSetSegmentationSurfaceInfo>(new LevelSetSegmentationSurfaceInfo(surface, featureAngle, surfaceSmoothing));
2500 }
2501}
2502
2507const SurfaceSkdTree & LevelSetSegmentationObject::getSearchTree() const {
2508 return m_surfaceInfo->getSearchTree();
2509}
2510
2515double LevelSetSegmentationObject::getFeatureAngle() const {
2516 return m_surfaceInfo->getFeatureAngle();
2517}
2518
2524bitpit::LevelSetSurfaceSmoothing LevelSetSegmentationObject::getSurfaceSmoothing() const {
2525 return m_surfaceInfo->getSurfaceSmoothing();
2526}
2527
2538void LevelSetSegmentationObject::fillCellLocationCache()
2539{
2540 // Cartesian patches are handled separately
2541 if (dynamic_cast<LevelSetCartesianKernel*>(m_kernel)) {
2542 fillCartesianCellZoneCache();
2543 return;
2544 }
2545
2546 // All other patches are handled with the base method.
2547 LevelSetObject::fillCellLocationCache();
2548}
2549
2562void LevelSetSegmentationObject::fillCellLocationCache(const std::vector<adaption::Info> &adaptionData)
2563{
2564 // Cartesian patches are handled separately
2565 //
2566 // Update is not implemented for Cartesian patches, the cells that should be inserted in
2567 // the cache will be evaluated considering all the mash not just the elements newly added.
2568 if (dynamic_cast<LevelSetCartesianKernel*>(m_kernel)) {
2569 fillCartesianCellZoneCache();
2570 return;
2571 }
2572
2573 // All other patches are handled with the base method
2574 LevelSetObject::fillCellLocationCache(adaptionData);
2575}
2576
2588void LevelSetSegmentationObject::fillCartesianCellZoneCache()
2589{
2590 // The function needs a Cartesian kernel
2591 const LevelSetCartesianKernel *cartesianKernel = dynamic_cast<LevelSetCartesianKernel*>(m_kernel);
2592 if (!dynamic_cast<LevelSetCartesianKernel*>(m_kernel)) {
2593 throw std::runtime_error("The function needs a Cartesian kernels.");
2594 }
2595
2596 // Get mesh information
2597 const VolCartesian &mesh = *(cartesianKernel->getMesh() ) ;
2598 int meshDimension = mesh.getDimension();
2599 long nCells = mesh.getCellCount();
2600
2601 std::array<double,3> meshMinPoint;
2602 std::array<double,3> meshMaxPoint;
2603 mesh.getBoundingBox(meshMinPoint, meshMaxPoint) ;
2604
2605 // Get surface information
2606 const SurfUnstructured &surface = getSurface();
2607
2608 // Initialize process list
2609 //
2610 // Process list is initialized with cells that are certainly inside the
2611 // narrow band. Those cells are the ones that contain the vertices of the
2612 // segments or the intersection between the segments and the bounding box
2613 // of the patch.
2614 std::unordered_set<long> processList;
2615
2616 std::vector<std::array<double,3>> intersectionPoints;
2617 std::vector<std::array<double,3>> segmentVertexCoords;
2618 for (const Cell &segment : surface.getCells()) {
2619 // Get segment info
2620 //
2621 // Since vertex information will be passed to the CG module, we need to get the
2622 // vertices in the proper order.
2623 ConstProxyVector<long> segmentVertexIds = surface.getFacetOrderedVertexIds(segment);
2624 std::size_t nSegmentVertices = segmentVertexIds.size();
2625
2626 // Get segment coordinates
2627 segmentVertexCoords.resize(nSegmentVertices);
2628 surface.getVertexCoords(nSegmentVertices, segmentVertexIds.data(), segmentVertexCoords.data());
2629
2630 // Add to the process list the cells that contain the vertices of the
2631 // segment or the intersection between the segment and the bounding box
2632 // of the patch.
2633 int nInnerVertices = 0;
2634 for (const std::array<double,3> &vertexPoint : segmentVertexCoords) {
2635 long cellId = mesh.locatePoint(vertexPoint);
2636 if (cellId < 0) {
2637 continue;
2638 }
2639
2640 processList.insert(cellId);
2641 ++nInnerVertices;
2642 }
2643
2644 if (nInnerVertices == 0) {
2645 if (CGElem::intersectBoxPolygon(meshMinPoint, meshMaxPoint, segmentVertexCoords, false, true, true, intersectionPoints, meshDimension)) {
2646 for (const std::array<double,3> &intersectionPoint : intersectionPoints){
2647 long cellId = mesh.locateClosestCell(intersectionPoint);
2648 processList.insert(cellId);
2649 }
2650 }
2651 }
2652 }
2653
2654 // Get cache for zone identification
2655 CellCacheCollection::ValueCache<char> *locationCache = getCellCache<char>(m_cellLocationCacheId);
2656
2657 // Cells with an unknown region are in the bulk
2658 for (long cellId = 0; cellId < nCells; ++cellId) {
2659 locationCache->insertEntry(cellId, static_cast<char>(LevelSetCellLocation::UNKNOWN));
2660 }
2661
2662 // Start filling the list of cells within the narrow band
2663 //
2664 // The initial process list is gradually expanded considering all the neighbours
2665 // inside the narrow band.
2666 std::vector<long> intersectedCellIds;
2667 std::vector<bool> alreadyProcessed(nCells, false);
2668 while (!processList.empty()) {
2669 // Get the cell to process
2670 long cellId = *(processList.begin());
2671 processList.erase(processList.begin());
2672
2673 // Skip cells already processed
2674 if (alreadyProcessed[cellId]) {
2675 continue;
2676 }
2677 alreadyProcessed[cellId] = true;
2678
2679 // Fill location cache for cells geometrically inside the narrow band
2680 LevelSetCellLocation cellLocation = fillCellGeometricNarrowBandLocationCache(cellId);
2681
2682 // Track intersected cells
2683 if (cellLocation == LevelSetCellLocation::NARROW_BAND_INTERSECTED) {
2684 intersectedCellIds.push_back(cellId);
2685 }
2686
2687 // Add unprocessed neighbours of cells inside the narrow band to the process list
2688 if (cellLocation != LevelSetCellLocation::UNKNOWN) {
2689 auto neighProcessor = [&processList, &alreadyProcessed](long neighId, int layer) {
2690 BITPIT_UNUSED(layer);
2691
2692 // Skip neighbours already processed
2693 if (alreadyProcessed[neighId]) {
2694 return false;
2695 }
2696
2697 // Update the process list
2698 processList.insert(neighId);
2699
2700 // Continue processing the other neighbours
2701 return false;
2702 };
2703
2704 mesh.processCellFaceNeighbours(cellId, 1, neighProcessor);
2705 }
2706 }
2707
2708 // Identify neighbours of cells that intersect the surface
2709 for (std::size_t cellId : intersectedCellIds) {
2710 // Process face neighbours
2711 auto neighProcessor = [this, &locationCache](long neighId, int layer) {
2712 BITPIT_UNUSED(layer);
2713
2714 // Skip neighbours whose region has already been identified
2715 if (getCellLocation(neighId) != LevelSetCellLocation::UNKNOWN) {
2716 return false;
2717 }
2718
2719 // The neighbour is inside the narrow band
2720 locationCache->insertEntry(neighId, static_cast<char>(LevelSetCellLocation::NARROW_BAND_NEIGHBOUR));
2721
2722 // Continue processing the other neighbours
2723 return false;
2724 };
2725
2726 mesh.processCellFaceNeighbours(cellId, 1, neighProcessor);
2727 }
2728
2729 // Cells with an unknown region are in the bulk
2730 for (long cellId = 0; cellId < nCells; ++cellId) {
2731 CellCacheCollection::ValueCache<char>::Entry locationCacheEntry = locationCache->findEntry(cellId);
2732 if (*locationCacheEntry == static_cast<char>(LevelSetCellLocation::UNKNOWN)) {
2733 locationCache->insertEntry(cellId, static_cast<char>(LevelSetCellLocation::BULK));
2734 }
2735 }
2736
2737#if BITPIT_ENABLE_MPI==1
2738 // Exchange ghost data
2739 if (mesh.isPartitioned()) {
2740 std::unique_ptr<DataCommunicator> dataCommunicator = m_kernel->createDataCommunicator();
2741 startCellCacheExchange(mesh.getGhostCellExchangeSources(), m_cellLocationCacheId, dataCommunicator.get());
2742 completeCellCacheExchange(mesh.getGhostCellExchangeTargets(), m_cellLocationCacheId, dataCommunicator.get());
2743 }
2744#endif
2745}
2746
2761LevelSetCellLocation LevelSetSegmentationObject::fillCellGeometricNarrowBandLocationCache(long id)
2762{
2763 // Early return if the cell is geometrically outside the narrow band
2764 double searchRadius = std::max(m_kernel->computeCellBoundingRadius(id), m_narrowBandSize);
2765 long cellSupport = evalCellSupport(id, searchRadius);
2766 if (cellSupport < 0) {
2767 return LevelSetCellLocation::UNKNOWN;
2768 }
2769
2770 // Evaluate levelset value
2771 std::array<double,3> cellCentroid = m_kernel->computeCellCentroid(id);
2772
2773 double cellCacheValue = _evalValue(cellCentroid, cellSupport, CELL_CACHE_IS_SIGNED);
2774 double cellUnsigendValue = std::abs(cellCacheValue);
2775
2776 // Update the cell location cache
2777 //
2778 // First we need to check if the cell intersects the surface, and only if it
2779 // doesn't we should check if its distance is lower than the narrow band size.
2780 //
2781 // When high order smoothing is active, there may be cells whose support is within the search
2782 // radius, but they are not intersected and their distance is less than the narrow band size.
2783 // These cells are not geometrically inside the narrow band, they are neighbours of cells
2784 // geometrically inside the narrow band and as such it's up to the caller of this function to
2785 // identify their cell location.
2786 LevelSetCellLocation cellLocation = LevelSetCellLocation::UNKNOWN;
2787 if (_intersectSurface(id, cellUnsigendValue, CELL_LOCATION_INTERSECTION_MODE) == LevelSetIntersectionStatus::TRUE) {
2788 cellLocation = LevelSetCellLocation::NARROW_BAND_INTERSECTED;
2789 } else if (cellUnsigendValue <= m_narrowBandSize) {
2790 cellLocation = LevelSetCellLocation::NARROW_BAND_DISTANCE;
2791 }
2792 assert((getSurfaceSmoothing() == LevelSetSurfaceSmoothing::HIGH_ORDER) || (cellLocation != LevelSetCellLocation::UNKNOWN));
2793
2794 if (cellLocation != LevelSetCellLocation::UNKNOWN) {
2795 CellCacheCollection::ValueCache<char> *locationCache = getCellCache<char>(m_cellLocationCacheId);
2796 locationCache->insertEntry(id, static_cast<char>(cellLocation));
2797 }
2798
2799 // Fill the cache of the evaluated fields
2800 //
2801 // Now that the cell location has been identified, we can fill the cache of the
2802 // evaluated fields.
2803 fillFieldCellCache(LevelSetField::SUPPORT, id, cellSupport);
2804 fillFieldCellCache(LevelSetField::VALUE, id, cellCacheValue);
2805
2806 // Return the location
2807 return cellLocation;
2808}
2809
2816const SurfUnstructured & LevelSetSegmentationObject::_evalCellSurface(long id) const
2817{
2818 BITPIT_UNUSED(id);
2819
2820 return getSurface();
2821}
2822
2829short LevelSetSegmentationObject::_evalCellSign(long id) const
2830{
2831 long support = evalCellSupport(id);
2832 std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
2833
2834 return _evalSign(centroid, support);
2835}
2836
2844double LevelSetSegmentationObject::_evalCellValue(long id, bool signedLevelSet) const
2845{
2846 long support = evalCellSupport(id);
2847 std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
2848
2849 return _evalValue(centroid, support, signedLevelSet);
2850}
2851
2859std::array<double,3> LevelSetSegmentationObject::_evalCellGradient(long id, bool signedLevelSet) const
2860{
2861 long support = evalCellSupport(id);
2862 std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
2863
2864 return _evalGradient(centroid, support, signedLevelSet);
2865}
2866
2879long LevelSetSegmentationObject::_evalCellSupport(long id, double searchRadius) const
2880{
2881 std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
2882
2883 return _evalSupport(centroid, searchRadius);
2884}
2885
2893std::array<double,3> LevelSetSegmentationObject::_evalCellNormal(long id, bool signedLevelSet) const
2894{
2895 long support = evalCellSupport(id);
2896 std::array<double,3> centroid = m_kernel->computeCellCentroid(id);
2897
2898 std::array<double, 3> projectionPoint;
2899 std::array<double, 3> projectionNormal;
2900 _evalProjection(centroid, support, signedLevelSet, &projectionPoint, &projectionNormal);
2901
2902 return projectionNormal;
2903}
2904
2911const SurfUnstructured & LevelSetSegmentationObject::_evalSurface(const std::array<double,3> &point) const
2912{
2913 BITPIT_UNUSED(point);
2914
2915 return getSurface();
2916}
2917
2924short LevelSetSegmentationObject::_evalSign(const std::array<double,3> &point) const
2925{
2926 long support = evalSupport(point);
2927
2928 return _evalSign(point, support);
2929}
2930
2938double LevelSetSegmentationObject::_evalValue(const std::array<double,3> &point, bool signedLevelSet) const
2939{
2940 long support = evalSupport(point);
2941
2942 return _evalValue(point, support, signedLevelSet);
2943}
2944
2952std::array<double,3> LevelSetSegmentationObject::_evalGradient(const std::array<double,3> &point, bool signedLevelSet) const
2953{
2954 long support = evalSupport(point);
2955
2956 return _evalGradient(point, support, signedLevelSet);
2957}
2958
2970void LevelSetSegmentationObject::_evalProjection(const std::array<double,3> &point,
2971 bool signedLevelSet,
2972 std::array<double, 3> *projectionPoint,
2973 std::array<double, 3> *projectionNormal) const
2974{
2975 // Get closest segment
2976 long support = evalSupport(point);
2977
2978 return _evalProjection(point, support, signedLevelSet, projectionPoint, projectionNormal);
2979}
2980
2988short LevelSetSegmentationObject::_evalSign(const std::array<double,3> &point, long support) const
2989{
2990 // Throw an error if the support is not valid
2991 if (support < 0) {
2992 if (empty()) {
2993 return levelSetDefaults::SIGN;
2994 }
2995
2996 throw std::runtime_error("Unable to evaluate the sign: the support is not valid.");
2997 }
2998
2999 LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support);
3000
3001 return evalValueSign(m_surfaceInfo->evalDistance(point, supportItr, true));
3002}
3003
3012double LevelSetSegmentationObject::_evalValue(const std::array<double,3> &point, long support,
3013 bool signedLevelSet) const
3014{
3015 // Early return if the support is not valid
3016 //
3017 // With an invalid support, only the unsigend levelset can be evaluated.
3018 if (support < 0) {
3019 if (!signedLevelSet || empty()) {
3020 return levelSetDefaults::VALUE;
3021 }
3022
3023 throw std::runtime_error("With an invalid support, only the unsigend levelset can be evaluated.");
3024 }
3025
3026 // Evaluate the distance of the point from the surface
3027 LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support);
3028 double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet);
3029
3030 // Early return if the point lies on the surface
3031 if (evalValueSign(distance) == 0) {
3032 return 0.;
3033 }
3034
3035 // Evaluate levelset value
3036 double value;
3037 if (signedLevelSet) {
3038 value = distance;
3039 } else {
3040 value = std::abs(distance);
3041 }
3042
3043 return value;
3044}
3045
3054std::array<double,3> LevelSetSegmentationObject::_evalGradient(const std::array<double,3> &point, long support,
3055 bool signedLevelSet) const
3056{
3057 // Early return if the support is not valid
3058 //
3059 // With an invalid support, only the unsigend levelset can be evaluated.
3060 if (support < 0) {
3061 if (!signedLevelSet || empty()) {
3062 return levelSetDefaults::GRADIENT;
3063 }
3064
3065 throw std::runtime_error("With an invalid support, only the unsigend levelset can be evaluated.");
3066 }
3067
3068 // Evaluate the distance of the point from the surface
3069 LevelSetSegmentationSurfaceInfo::SegmentConstIterator supportItr = getSurface().getCellConstIterator(support);
3070 std::array<double,3> distanceVector;
3071 double distance = m_surfaceInfo->evalDistance(point, supportItr, signedLevelSet, &distanceVector);
3072
3073 // Early return if the point lies on the surface
3074 if (evalValueSign(distance) == 0) {
3075 if (signedLevelSet) {
3076 return m_surfaceInfo->evalNormal(point, supportItr);
3077 } else {
3078 return {{0., 0., 0.}};
3079 }
3080 }
3081
3082 // Evaluate levelset gradient
3083 std::array<double,3> gradient = distanceVector / distance;
3084
3085 return gradient;
3086}
3087
3094long LevelSetSegmentationObject::_evalSupport(const std::array<double,3> &point) const
3095{
3096 return _evalSupport(point, std::numeric_limits<double>::max());
3097}
3098
3107long LevelSetSegmentationObject::_evalSupport(const std::array<double,3> &point, double searchRadius) const
3108{
3109 long closestSegmentId;
3110 double closestDistance;
3111 m_surfaceInfo->getSearchTree().findPointClosestCell(point, searchRadius, &closestSegmentId, &closestDistance);
3112
3113 return closestSegmentId;
3114}
3115
3128void LevelSetSegmentationObject::_evalProjection(const std::array<double,3> &point,
3129 long support,
3130 bool signedLevelSet,
3131 std::array<double, 3> *projectionPoint,
3132 std::array<double, 3> *projectionNormal) const
3133{
3134 // Early return if the support is not valid
3135 //
3136 // With an invalid support, only the unsigend levelset can be evaluated.
3137 if (support < 0) {
3138 if (!signedLevelSet || empty()) {
3139 (*projectionPoint) = levelSetDefaults::POINT;
3140 (*projectionNormal) = levelSetDefaults::NORMAL;
3141 }
3142
3143 throw std::runtime_error("With an invalid support, only the unsigend levelset can be evaluated.");
3144 }
3145
3146 // Get closest segment
3147 LevelSetSegmentationSurfaceInfo::SegmentConstIterator segmentItr = getSurface().getCellConstIterator(support);
3148
3149 // Eval projection point and normal
3150 m_surfaceInfo->evalProjection(point, segmentItr, projectionPoint, projectionNormal);
3151
3152 // If an unsigned evaluation is requested, the orientation of the surface should be discarded
3153 // and in order to have a normal that is agnostic with respect the two sides of the surface.
3154 if (!signedLevelSet) {
3155 (*projectionNormal) *= static_cast<double>(evalSign(point));
3156 }
3157}
3158
3165double LevelSetSegmentationObject::getMinSurfaceFeatureSize() const {
3166
3167 const SurfUnstructured &surface = getSurface();
3168
3169 bool minimumValid = false;
3170 double minimumSize = levelSetDefaults::SIZE;
3171 for (const Cell &cell : surface.getCells()) {
3172 double segmentSize = surface.evalCellSize(cell.getId());
3173 if (segmentSize < 0) {
3174 continue;
3175 }
3176
3177 minimumValid = true;
3178 minimumSize = std::min(segmentSize, minimumSize);
3179 }
3180
3181 if (!minimumValid) {
3182 minimumSize = - levelSetDefaults::SIZE;
3183 }
3184
3185 return minimumSize;
3186}
3187
3194double LevelSetSegmentationObject::getMaxSurfaceFeatureSize() const {
3195
3196 const SurfUnstructured &surface = getSurface();
3197
3198 double maximumSize = - levelSetDefaults::SIZE;
3199 for (const Cell &cell : surface.getCells()) {
3200 double segmentSize = surface.evalCellSize(cell.getId());
3201 maximumSize = std::max(segmentSize, maximumSize);
3202 }
3203
3204 return maximumSize;
3205}
3206
3217
3225LevelSetBooleanObject<LevelSetSegmentationBaseObject>::LevelSetBooleanObject( int id, LevelSetBooleanOperation op, const std::vector<const LevelSetSegmentationBaseObject *> &sourceObjects )
3226 : LevelSetBooleanBaseObject<LevelSetSegmentationBaseObject>(id, op, sourceObjects) {
3227}
3228
3236
3243const SurfUnstructured & LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellSurface(long id) const
3244{
3245 return getCellReferenceObject(id)->evalCellSurface(id);
3246}
3247
3254long LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellSupport(long id, double searchRadius) const
3255{
3256 return getCellReferenceObject(id)->evalCellSupport(id, searchRadius);
3257}
3258
3265int LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellPart(long id) const
3266{
3267 return getCellReferenceObject(id)->evalCellPart(id);
3268}
3269
3277std::array<double,3> LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalCellNormal(long id, bool signedLevelSet) const
3278{
3279 return _evalCellFunction<std::array<double,3>>(id, signedLevelSet, [&id, signedLevelSet] (const LevelSetBooleanResult<LevelSetSegmentationBaseObject> &result)
3280 {
3281 const LevelSetSegmentationBaseObject *resultObject = result.getObject();
3282 if ( !resultObject ) {
3283 return levelSetDefaults::NORMAL;
3284 }
3285
3286 std::array<double,3> normal = resultObject->evalCellNormal(id, signedLevelSet);
3287 if (signedLevelSet) {
3288 normal *= static_cast<double>(result.getObjectSign());
3289 }
3290
3291 return normal;
3292 });
3293}
3294
3301const SurfUnstructured & LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalSurface(const std::array<double,3> &point) const
3302{
3303 return getReferenceObject(point)->evalSurface(point);
3304}
3305
3312long LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point) const
3313{
3314 return getReferenceObject(point)->evalSupport(point);
3315}
3316
3325long LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point, double searchRadius) const
3326{
3327 return getReferenceObject(point)->evalSupport(point, searchRadius);
3328}
3329
3341void LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalProjection(const std::array<double,3> &point,
3342 bool signedLevelSet,
3343 std::array<double, 3> *projectionPoint,
3344 std::array<double, 3> *projectionNormal) const
3345{
3346 return _evalFunction<void>(point, signedLevelSet, [&point, projectionPoint, projectionNormal, signedLevelSet] (const LevelSetBooleanResult<LevelSetSegmentationBaseObject> &result)
3347 {
3348 const LevelSetSegmentationBaseObject *resultObject = result.getObject();
3349 if ( !resultObject ) {
3350 (*projectionNormal) = levelSetDefaults::NORMAL;
3351 (*projectionPoint) = levelSetDefaults::POINT;
3352 }
3353
3354 resultObject->evalProjection(point, signedLevelSet, projectionPoint, projectionNormal);
3355 if (signedLevelSet) {
3356 (*projectionNormal) *= static_cast<double>(result.getObjectSign());
3357 }
3358 });
3359}
3360
3367int LevelSetBooleanObject<LevelSetSegmentationBaseObject>::_evalPart(const std::array<double,3> &point) const
3368{
3369 return getReferenceObject(point)->evalPart(point);
3370}
3371
3382
3390
3397const SurfUnstructured & LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellSurface(long id) const
3398{
3399 return getSourceObject()->evalCellSurface(id);
3400}
3401
3414long LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellSupport(long id, double searchRadius) const
3415{
3416 return getSourceObject()->evalCellSupport(id, searchRadius);
3417}
3418
3425int LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellPart(long id) const
3426{
3427 return getCellReferenceObject(id)->evalCellPart(id);
3428}
3429
3437std::array<double,3> LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalCellNormal(long id, bool signedLevelSet) const
3438{
3439 std::array<double,3> normal = getSourceObject()->evalCellNormal(id, signedLevelSet);
3440 if (signedLevelSet) {
3441 normal *= -1.;
3442 }
3443
3444 return normal;
3445}
3446
3453const SurfUnstructured & LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalSurface(const std::array<double,3> &point) const
3454{
3455 return getSourceObject()->evalSurface(point);
3456}
3457
3464long LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point) const
3465{
3466 return getSourceObject()->evalSupport(point);
3467}
3468
3477long LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalSupport(const std::array<double,3> &point, double searchRadius) const
3478{
3479 return getSourceObject()->evalSupport(point, searchRadius);
3480}
3481
3488int LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalPart(const std::array<double,3> &point) const
3489{
3490 return getSourceObject()->evalPart(point);
3491}
3492
3504void LevelSetComplementObject<LevelSetSegmentationBaseObject>::_evalProjection(const std::array<double,3> &point,
3505 bool signedLevelSet,
3506 std::array<double, 3> *projectionPoint,
3507 std::array<double, 3> *projectionNormal) const
3508{
3509 getSourceObject()->evalProjection(point, signedLevelSet, projectionPoint, projectionNormal);
3510 if (signedLevelSet) {
3511 (*projectionNormal) *= -1.;
3512 }
3513}
3514
3515}
bitpit::Cell
The Cell class defines the cells.
Definition cell.hpp:42
bitpit::Element::getId
long getId() const
Definition element.cpp:510
bitpit::Element::getType
ElementType getType() const
Definition element.cpp:551
bitpit::Element::getVertexIds
static ConstProxyVector< long > getVertexIds(ElementType type, const long *connectivity)
Definition element.cpp:1193
bitpit::Element::getVertexCount
int getVertexCount() const
Definition element.cpp:1152
bitpit::Element::getPID
int getPID() const
Definition element.cpp:581
bitpit::LevelSetBooleanBaseObject
Base class which deals with boolean operation between two LevelSetObjects.
Definition levelSetBooleanObject.hpp:68
bitpit::LevelSetBooleanObject
Class which deals with boolean operation between two LevelSetObjects.
Definition levelSetBooleanObject.hpp:113
bitpit::LevelSetBooleanResult
Allow to evaluate the result of a boolean operation between two LevelSetObjects.
Definition levelSetBooleanObject.hpp:44
bitpit::LevelSetBooleanResult::getObject
const SourceLevelSetObject * getObject() const
Definition levelSetBooleanObject.tpp:110
bitpit::LevelSetCartesianKernel
Implements LevelSetKernel for cartesian meshes.
Definition levelSetCartesianKernel.hpp:34
bitpit::LevelSetCartesianKernel::getMesh
VolCartesian * getMesh() const override
Definition levelSetCartesianKernel.cpp:69
bitpit::LevelSetComplementBaseObject
Base class that allows to evaluate the complement of a LevelSetObjects.
Definition levelSetComplementObject.hpp:38
bitpit::LevelSetComplementObject
Class that allows to evaluate the complement of a LevelSetObjects.
Definition levelSetComplementObject.hpp:71
bitpit::LevelSetKernel::getMesh
virtual VolumeKernel * getMesh() const
Definition levelSetKernel.cpp:83
bitpit::LevelSetKernel::createDataCommunicator
std::unique_ptr< DataCommunicator > createDataCommunicator() const
Definition levelSetKernel.cpp:274
bitpit::LevelSetObject::CELL_LOCATION_INTERSECTION_MODE
static const LevelSetIntersectionMode CELL_LOCATION_INTERSECTION_MODE
Definition levelSetObject.hpp:132
bitpit::LevelSetObject::fillFieldCellCache
void fillFieldCellCache(LevelSetField field, const std::vector< long > &cellIds)
Definition levelSetObject.cpp:3043
bitpit::LevelSetObject::evalValueSign
short evalValueSign(double value) const
Definition levelSetObject.cpp:1669
bitpit::LevelSetObject::m_kernel
LevelSetKernel * m_kernel
Definition levelSetObject.hpp:134
bitpit::LevelSetObject::completeCellCacheExchange
void completeCellCacheExchange(const std::unordered_map< int, std::vector< long > > &sendCellIds, std::size_t cacheIds, DataCommunicator *)
Definition levelSetObject.cpp:2190
bitpit::LevelSetObject::addVTKOutputData
virtual void addVTKOutputData(LevelSetField field, const std::string &objectName)
Definition levelSetObject.cpp:1990
bitpit::LevelSetObject::CELL_CACHE_IS_SIGNED
static const bool CELL_CACHE_IS_SIGNED
Definition levelSetObject.hpp:131
bitpit::LevelSetObject::m_narrowBandSize
double m_narrowBandSize
Size of narrow band.
Definition levelSetObject.hpp:140
bitpit::LevelSetObject::getVTKOutputFieldName
virtual std::string getVTKOutputFieldName(LevelSetField field) const
Definition levelSetObject.cpp:2040
bitpit::LevelSetObject::evalSign
virtual short evalSign(const std::array< double, 3 > &point) const
Definition levelSetObject.cpp:1621
bitpit::LevelSetObject::_intersectSurface
virtual LevelSetIntersectionStatus _intersectSurface(long, double distance, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const
Definition levelSetObject.cpp:1402
bitpit::LevelSetObject::getSupportedFields
virtual LevelSetFieldset getSupportedFields() const
Definition levelSetObject.cpp:175
bitpit::LevelSetObject::fillCellLocationCache
virtual void fillCellLocationCache()
Definition levelSetObject.cpp:550
bitpit::LevelSetObject::evalCellSign
virtual short evalCellSign(long id) const
Definition levelSetObject.cpp:1481
bitpit::LevelSetObject::getCellLocation
LevelSetCellLocation getCellLocation(long id) const
Definition levelSetObject.cpp:483
bitpit::LevelSetObject::m_cellLocationCacheId
std::size_t m_cellLocationCacheId
Id of the cache that will keep track if cell zones.
Definition levelSetObject.hpp:142
bitpit::LevelSetObject::flushVTKOutputData
virtual void flushVTKOutputData(std::fstream &stream, VTKFormat format, LevelSetField field) const
Definition levelSetObject.cpp:2095
bitpit::LevelSetObject::createFieldCellCache
virtual std::size_t createFieldCellCache(LevelSetField field, std::size_t cacheId=CellCacheCollection::NULL_CACHE_ID)
Definition levelSetObject.cpp:2842
bitpit::LevelSetObject::getVTKOutputDataName
std::string getVTKOutputDataName(LevelSetField field, const std::string &objectName) const
Definition levelSetObject.cpp:2024
bitpit::LevelSetObject::startCellCacheExchange
void startCellCacheExchange(const std::unordered_map< int, std::vector< long > > &recvCellIds, std::size_t cacheIds, DataCommunicator *) const
Definition levelSetObject.cpp:2139
bitpit::LevelSetSegmentationBaseObject
Implements visitor pattern fo segmentated geometries.
Definition levelSetSegmentationObject.hpp:126
bitpit::LevelSetSegmentationBaseObject::evalProjection
void evalProjection(const std::array< double, 3 > &point, bool signedLevelSet, std::array< double, 3 > *projectionPoint, std::array< double, 3 > *projectionNormal) const
Definition levelSetSegmentationObject.cpp:2037
bitpit::LevelSetSegmentationBaseObject::fillFieldCellCache
void fillFieldCellCache(LevelSetField field, long id) override
Definition levelSetSegmentationObject.cpp:1720
bitpit::LevelSetSegmentationBaseObject::flushVTKOutputData
void flushVTKOutputData(std::fstream &stream, VTKFormat format, LevelSetField field) const override
Definition levelSetSegmentationObject.cpp:2141
bitpit::LevelSetSegmentationBaseObject::evalNormal
std::array< double, 3 > evalNormal(const std::array< double, 3 > &point, bool signedLevelSet) const
Definition levelSetSegmentationObject.cpp:1992
bitpit::LevelSetSegmentationBaseObject::evalSurface
const SurfUnstructured & evalSurface(const std::array< double, 3 > &point) const
Definition levelSetSegmentationObject.cpp:1969
bitpit::LevelSetSegmentationBaseObject::getSupportedFields
LevelSetFieldset getSupportedFields() const override
Definition levelSetSegmentationObject.cpp:1704
bitpit::LevelSetSegmentationBaseObject::getNormal
std::array< double BITPIT_COMMA 3 > getNormal(long cellId) const
Definition levelSetSegmentationObject.cpp:2194
bitpit::LevelSetSegmentationBaseObject::_intersectSurface
LevelSetIntersectionStatus _intersectSurface(long, double distance, LevelSetIntersectionMode=LevelSetIntersectionMode::FAST_FUZZY) const override
Definition levelSetSegmentationObject.cpp:1806
bitpit::LevelSetSegmentationBaseObject::_evalPart
virtual int _evalPart(const std::array< double, 3 > &point) const
Definition levelSetSegmentationObject.cpp:2066
bitpit::LevelSetSegmentationBaseObject::addVTKOutputData
void addVTKOutputData(LevelSetField field, const std::string &objectName) override
Definition levelSetSegmentationObject.cpp:2081
bitpit::LevelSetSegmentationBaseObject::evalCellNormal
std::array< double, 3 > evalCellNormal(long id, bool signedLevelSet) const
Definition levelSetSegmentationObject.cpp:1827
bitpit::LevelSetSegmentationBaseObject::createFieldCellCache
std::size_t createFieldCellCache(LevelSetField field, std::size_t cacheId=CellCacheCollection::NULL_CACHE_ID) override
Definition levelSetSegmentationObject.cpp:1687
bitpit::LevelSetSegmentationBaseObject::evalPart
int evalPart(const std::array< double, 3 > &point) const
Definition levelSetSegmentationObject.cpp:1980
bitpit::LevelSetSegmentationBaseObject::evalCellSupport
long evalCellSupport(long id, double searchRadius=AUTOMATIC_SEARCH_RADIUS) const
Definition levelSetSegmentationObject.cpp:1875
bitpit::LevelSetSegmentationBaseObject::AUTOMATIC_SEARCH_RADIUS
static BITPIT_PUBLIC_API const double AUTOMATIC_SEARCH_RADIUS
Definition levelSetSegmentationObject.hpp:129
bitpit::LevelSetSegmentationBaseObject::evalSupport
long evalSupport(const std::array< double, 3 > &point) const
Definition levelSetSegmentationObject.cpp:2008
bitpit::LevelSetSegmentationBaseObject::evalCellSurface
const SurfUnstructured & evalCellSurface(long id) const
Definition levelSetSegmentationObject.cpp:1740
bitpit::LevelSetSegmentationBaseObject::getVTKOutputFieldName
std::string getVTKOutputFieldName(LevelSetField field) const override
Definition levelSetSegmentationObject.cpp:2113
bitpit::LevelSetSegmentationObject
Implements visitor pattern fo segmentated geometries.
Definition levelSetSegmentationObject.hpp:186
bitpit::LevelSetSegmentationObject::_evalSurface
const SurfUnstructured & _evalSurface(const std::array< double, 3 > &point) const override
Definition levelSetSegmentationObject.cpp:2911
bitpit::LevelSetSegmentationObject::clone
LevelSetSegmentationObject * clone() const override
Definition levelSetSegmentationObject.cpp:2361
bitpit::LevelSetSegmentationObject::_evalCellSupport
long _evalCellSupport(long id, double searchRadius=AUTOMATIC_SEARCH_RADIUS) const override
Definition levelSetSegmentationObject.cpp:2879
bitpit::LevelSetSegmentationObject::fillCellGeometricNarrowBandLocationCache
LevelSetCellLocation fillCellGeometricNarrowBandLocationCache(long id) override
Definition levelSetSegmentationObject.cpp:2761
bitpit::LevelSetSegmentationObject::_evalSupport
long _evalSupport(const std::array< double, 3 > &point) const override
Definition levelSetSegmentationObject.cpp:3094
bitpit::LevelSetSegmentationObject::_evalSign
short _evalSign(const std::array< double, 3 > &point) const override
Definition levelSetSegmentationObject.cpp:2924
bitpit::LevelSetSegmentationObject::_evalGradient
std::array< double, 3 > _evalGradient(const std::array< double, 3 > &point, bool signedLevelSet) const override
Definition levelSetSegmentationObject.cpp:2952
bitpit::LevelSetSegmentationObject::_evalCellValue
double _evalCellValue(long id, bool signedLevelSet) const override
Definition levelSetSegmentationObject.cpp:2844
bitpit::LevelSetSegmentationObject::_evalCellGradient
std::array< double, 3 > _evalCellGradient(long id, bool signedLevelSet) const override
Definition levelSetSegmentationObject.cpp:2859
bitpit::LevelSetSegmentationObject::getSurfaceSmoothing
LevelSetSurfaceSmoothing getSurfaceSmoothing() const
Definition levelSetSegmentationObject.cpp:2524
bitpit::LevelSetSegmentationObject::_evalValue
double _evalValue(const std::array< double, 3 > &point, bool signedLevelSet) const override
Definition levelSetSegmentationObject.cpp:2938
bitpit::LevelSetSegmentationObject::setSurface
void setSurface(std::unique_ptr< const SurfUnstructured > &&surface, bool force=false)
Definition levelSetSegmentationObject.cpp:2387
bitpit::LevelSetSegmentationObject::getSurface
const SurfUnstructured & getSurface() const
Definition levelSetSegmentationObject.cpp:2369
bitpit::LevelSetSegmentationObject::_evalCellNormal
std::array< double, 3 > _evalCellNormal(long id, bool signedLevelSet) const override
Definition levelSetSegmentationObject.cpp:2893
bitpit::LevelSetSegmentationObject::_evalCellSurface
const SurfUnstructured & _evalCellSurface(long id) const override
Definition levelSetSegmentationObject.cpp:2816
bitpit::LevelSetSegmentationObject::_evalProjection
void _evalProjection(const std::array< double, 3 > &point, bool signedLevelSet, std::array< double, 3 > *projectionPoint, std::array< double, 3 > *projectionNormal) const override
Definition levelSetSegmentationObject.cpp:2970
bitpit::LevelSetSegmentationSurfaceInfo::evalProjection
void evalProjection(const std::array< double, 3 > &point, const SegmentConstIterator &segmentItr, std::array< double, 3 > *projectionPoint, std::array< double, 3 > *projectionNormal) const
Definition levelSetSegmentationObject.cpp:1399
bitpit::LevelSetSegmentationSurfaceInfo::evalDistance
double evalDistance(const std::array< double, 3 > &point, const SegmentConstIterator &segmentItr, bool signedDistance, std::array< double, 3 > *distanceVector) const
Definition levelSetSegmentationObject.cpp:231
bitpit::LevelSetSegmentationSurfaceInfo::evalDistanceVector
std::array< double, 3 > evalDistanceVector(const std::array< double, 3 > &point, const SegmentConstIterator &segmentItr) const
Definition levelSetSegmentationObject.cpp:290
bitpit::LevelSetSegmentationSurfaceInfo::setSurface
void setSurface(std::unique_ptr< const SurfUnstructured > &&surface, double featureAngle=DEFAULT_FEATURE_ANGLE)
Definition levelSetSegmentationObject.cpp:134
bitpit::LevelSetSegmentationSurfaceInfo::evalNormal
std::array< double, 3 > evalNormal(const std::array< double, 3 > &point, const SegmentConstIterator &segmentItr) const
Definition levelSetSegmentationObject.cpp:307
bitpit::LevelSetSegmentationSurfaceInfo::DEFAULT_FEATURE_ANGLE
static BITPIT_PUBLIC_API const double DEFAULT_FEATURE_ANGLE
Definition levelSetSegmentationObject.hpp:57
bitpit::LevelSetSegmentationSurfaceInfo::getSurfaceSmoothing
LevelSetSurfaceSmoothing getSurfaceSmoothing() const
Definition levelSetSegmentationObject.cpp:218
bitpit::LevelSetSegmentationSurfaceInfo::evalPseudoNormal
std::array< double, 3 > evalPseudoNormal(const SurfUnstructured::CellConstIterator &segmentIterator, const double *lambda) const
Definition levelSetSegmentationObject.cpp:325
bitpit::LevelSetValueCache
The class LevelSetCache is the base class for defining caches that store values.
Definition levelSetCache.hpp:155
bitpit::PatchKernel::getAdjacenciesBuildStrategy
AdjacenciesBuildStrategy getAdjacenciesBuildStrategy() const
Definition patch_kernel.cpp:5785
bitpit::PatchKernel::getCells
PiercedVector< Cell > & getCells()
Definition patch_kernel.cpp:2551
bitpit::PatchKernel::empty
bool empty(bool global=true) const
Definition patch_kernel.cpp:2379
bitpit::PatchKernel::getVTK
VTKUnstructuredGrid & getVTK()
Definition patch_kernel.cpp:7784
bitpit::PatchKernel::getCell
Cell & getCell(long id)
Definition patch_kernel.cpp:2572
bitpit::PatchKernel::getVertices
PiercedVector< Vertex > & getVertices()
Definition patch_kernel.cpp:1664
bitpit::PatchKernel::processCellFaceNeighbours
void processCellFaceNeighbours(long seedId, int nLayers, Function function) const
Definition patch_kernel.tpp:497
bitpit::PatchKernel::findCellVertexNeighs
std::vector< long > findCellVertexNeighs(long id, bool complete=true) const
Definition patch_kernel.cpp:3778
bitpit::PatchKernel::getVertexCoords
const std::array< double, 3 > & getVertexCoords(long id) const
Definition patch_kernel.cpp:2338
bitpit::PatchKernel::getTol
double getTol() const
Definition patch_kernel.cpp:7596
bitpit::PatchKernel::getCellConstIterator
CellConstIterator getCellConstIterator(long id) const
Definition patch_kernel.cpp:2698
bitpit::PiercedStorageIterator
Iterator for the class PiercedStorage.
Definition piercedStorageIterator.hpp:53
bitpit::PiercedStorageSyncSlave::setStaticKernel
void setStaticKernel(const PiercedKernel< id_t > *kernel)
Definition piercedStorage.tpp:238
bitpit::PiercedStorageSyncSlave::unsetKernel
void unsetKernel(bool release=true)
Definition piercedStorage.tpp:310
bitpit::PiercedStorage::rawAt
__PS_REFERENCE__ rawAt(std::size_t pos, std::size_t offset=0)
Definition piercedStorage.tpp:1485
bitpit::PiercedStorage::end
iterator end() noexcept
Definition piercedStorage.tpp:1659
bitpit::PiercedStorage::rawData
__PS_POINTER__ rawData(std::size_t pos, std::size_t offset=0)
Definition piercedStorage.tpp:1239
bitpit::PiercedStorage::fill
void fill(const value_t &value)
Definition piercedStorage.tpp:1168
bitpit::PiercedStorage::begin
iterator begin() noexcept
Definition piercedStorage.tpp:1648
bitpit::ProxyVector
Metafunction for generating a list of elements that can be either stored in an external vectror or,...
Definition proxyVector.hpp:293
bitpit::SurfUnstructured
The SurfUnstructured class defines an unstructured surface triangulation.
Definition surfunstructured.hpp:36
bitpit::SurfaceKernel::getFacetOrderedLocalEdge
int getFacetOrderedLocalEdge(const Cell &facet, std::size_t n) const
Definition surface_kernel.cpp:2138
bitpit::SurfaceKernel::evalVertexNormals
virtual void evalVertexNormals(long id, int vertex, std::size_t nVertexNeighs, const long *vertexNeighs, double limit, std::array< double, 3 > *unlimitedNormal, std::array< double, 3 > *limitedNormal) const
Definition surface_kernel.cpp:893
bitpit::SurfaceKernel::evalFacetNormal
virtual std::array< double, 3 > evalFacetNormal(long, const std::array< double, 3 > &orientation={{0., 0., 1.}}) const
Definition surface_kernel.cpp:705
bitpit::SurfaceKernel::getFacetOrderedLocalVertex
int getFacetOrderedLocalVertex(const Cell &facet, std::size_t n) const
Definition surface_kernel.cpp:2052
bitpit::SurfaceKernel::getFacetOrderedVertexIds
ConstProxyVector< long > getFacetOrderedVertexIds(const Cell &facet) const
Definition surface_kernel.cpp:1992
bitpit::SurfaceKernel::evalLimitedVertexNormal
std::array< double, 3 > evalLimitedVertexNormal(long, int, double) const
Definition surface_kernel.cpp:834
bitpit::SurfaceKernel::evalBarycentricCoordinates
void evalBarycentricCoordinates(long id, const std::array< double, 3 > &point, double *lambda) const
Definition surface_kernel.cpp:263
bitpit::SurfaceKernel::evalEdgeNormals
virtual void evalEdgeNormals(long id, int edge, double limit, std::array< double, 3 > *unlimitedNormal, std::array< double, 3 > *limitedNormal) const
Definition surface_kernel.cpp:983
bitpit::SurfaceKernel::evalCellSize
double evalCellSize(long id) const override
Definition surface_kernel.cpp:232
bitpit::SurfaceSkdTree
The SurfaceSkdTree implements a Bounding Volume Hierarchy tree for surface patches.
Definition surface_skd_tree.hpp:33
bitpit::VTK
A base class for VTK input output.
Definition VTK.hpp:298
bitpit::VTK::addData
VTKField & addData(VTKField &&field)
Definition VTK.cpp:281
bitpit::VolumeKernel
The VolumeKernel class provides an interface for defining volume patches.
Definition volume_kernel.hpp:33
bitpit::CGElem::projectPointTriangle
array3D projectPointTriangle(array3D const &, array3D const &, array3D const &, array3D const &)
Definition CG_elem.cpp:1049
bitpit::CGElem::projectPointSegment
array3D projectPointSegment(array3D const &, array3D const &, array3D const &)
Definition CG_elem.cpp:996
bitpit::CGElem::convertBarycentricToFlagSegment
int convertBarycentricToFlagSegment(std::array< double, 2 > const &, double tolerance=DEFAULT_DISTANCE_TOLERANCE)
Definition CG_elem.cpp:618
bitpit::CGElem::convertBarycentricToFlagPolygon
int convertBarycentricToFlagPolygon(std::vector< double > const &, double tolerance=DEFAULT_DISTANCE_TOLERANCE)
Definition CG_elem.cpp:686
bitpit::CGElem::intersectBoxPolygon
bool intersectBoxPolygon(array3D const &, array3D const &, std::vector< array3D > const &, int dim=3, const double distanceTolerance=DEFAULT_DISTANCE_TOLERANCE)
Definition CG_elem.cpp:2284
bitpit::CGElem::projectPointPolygon
array3D projectPointPolygon(array3D const &, std::vector< array3D > const &)
Definition CG_elem.cpp:1214
sum
void sum(const std::array< T, d > &x, T1 &s)
Definition MathOperators_array.tpp:544
norm
double norm(const std::array< T, d > &x, int p)
Definition MathOperators_array.tpp:754
crossProduct
std::array< T, 3 > crossProduct(const std::array< T, 3 > &x, const std::array< T, 3 > &y)
Definition MathOperators_array.tpp:890
sign
T sign(const T &)
Definition Operators.tpp:38
dotProduct
T dotProduct(const std::array< T, d > &x, const std::array< T, d > &y)
Definition MathOperators_array.tpp:851
norm2
double norm2(const std::array< T, d > &x)
Definition MathOperators_array.tpp:714
bitpit::VTKFormat
VTKFormat
Definition VTK.hpp:92
BITPIT_UNUSED
#define BITPIT_UNUSED(variable)
Definition compiler.hpp:63
BITPIT_CREATE_WORKSPACE
#define BITPIT_CREATE_WORKSPACE(workspace, item_type, size, stack_size)
Definition commonUtils.hpp:76
bitpit::LevelSetIntersectionMode
LevelSetIntersectionMode
Definition levelSetCommon.hpp:118
bitpit::LevelSetField
LevelSetField
Definition levelSetCommon.hpp:167
bitpit::LevelSetBooleanOperation
LevelSetBooleanOperation
Definition levelSetCommon.hpp:96
bitpit::LevelSetIntersectionStatus
LevelSetIntersectionStatus
Definition levelSetCommon.hpp:107
bitpit::LevelSetSurfaceSmoothing
LevelSetSurfaceSmoothing
Definition levelSetCommon.hpp:184
bitpit::LevelSetCellLocation
LevelSetCellLocation
Definition levelSetCommon.hpp:80
bitpit::LevelSetField::NORMAL
@ NORMAL
bitpit::LevelSetField::SUPPORT
@ SUPPORT
bitpit::LevelSetField::PART
@ PART
bitpit::LevelSetField::VALUE
@ VALUE
bitpit::LevelSetIntersectionStatus::FALSE
@ FALSE
bitpit::LevelSetIntersectionStatus::TRUE
@ TRUE
bitpit::LevelSetSurfaceSmoothing::LOW_ORDER
@ LOW_ORDER
bitpit::LevelSetSurfaceSmoothing::HIGH_ORDER
@ HIGH_ORDER
bitpit::LevelSetCellLocation::UNKNOWN
@ UNKNOWN
Unknown location.
bitpit::LevelSetCellLocation::BULK
@ BULK
Bulk zone.
bitpit::LevelSetCellLocation::NARROW_BAND_NEIGHBOUR
@ NARROW_BAND_NEIGHBOUR
bitpit::LevelSetCellLocation::NARROW_BAND_INTERSECTED
@ NARROW_BAND_INTERSECTED
Narrow band zone, the cell intersects the surface.
bitpit::LevelSetCellLocation::NARROW_BAND_DISTANCE
@ NARROW_BAND_DISTANCE
bitpit::ElementType
ElementType
Definition element_type.hpp:30
bitpit::levelSetDefaults::GRADIENT
const std::array< double, 3 > GRADIENT
Definition levelSetCommon.hpp:48
bitpit::levelSetDefaults::POINT
const std::array< double, 3 > POINT
Definition levelSetCommon.hpp:49
bitpit::levelSetDefaults::NORMAL
const std::array< double, 3 > NORMAL
Definition levelSetCommon.hpp:55
--- layout: doxygen_footer ---
Generated on Wed Jun 19 2024 10:21:05 for bitpit by   doxygen 1.10.0