diff --git a/ugbase/lib_grid/CMakeLists.txt b/ugbase/lib_grid/CMakeLists.txt index eadce4b15..49a8af4d2 100644 --- a/ugbase/lib_grid/CMakeLists.txt +++ b/ugbase/lib_grid/CMakeLists.txt @@ -84,6 +84,8 @@ set(srcAlgorithms algorithms/debug_util.cpp algorithms/extrusion/cylinder_extrusion.cpp algorithms/extrusion/expand_layers.cpp algorithms/extrusion/expand_layers_arte.cpp + algorithms/extrusion/expand_layers_arte3D.cpp + algorithms/extrusion/ArteExpandFracs3D.cpp algorithms/projections/overlying_subset_finder.hpp algorithms/projections/z_ray_tracer.hpp refinement/projectors/projection_handler.cpp diff --git a/ugbase/lib_grid/algorithms/extrusion/ArteExpandFracs3D.cpp b/ugbase/lib_grid/algorithms/extrusion/ArteExpandFracs3D.cpp new file mode 100644 index 000000000..3d950716e --- /dev/null +++ b/ugbase/lib_grid/algorithms/extrusion/ArteExpandFracs3D.cpp @@ -0,0 +1,2799 @@ +/* + * ArteExpandFracs3D.cpp + * + * Created on: 06.10.2024 + * Author: Markus M. Knodel + * + * * expand fractures using the Arte algorithm, 3D case + * + * Author: Markus Knodel, inspired by Arte from Fuchs and Sebastian Reiters code for fracture expansion without Arte + * + * implementing a class that gives the basic tools for Arte in 3D + * might be templated at a later stage to fulfill 2D and 3D purposes, if suitable + * ( so far 2D case one entire function, not so perfect, but running....) + * + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + +#include + +#include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h" +#include "lib_grid/callbacks/callbacks.h" +#include "lib_grid/grid/grid_util.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "support.h" +#include "support3D.h" + + +#include + + + +namespace ug +{ + +ArteExpandFracs3D::ArteExpandFracs3D( + Grid & grid, SubsetHandler & sh, + std::vector const & fracInfos, + bool useTrianglesInDiamonds, bool establishDiamonds ) + : m_grid(grid), + m_sh(sh), + m_fracInfos(fracInfos), + m_useTrianglesInDiamonds(useTrianglesInDiamonds), + m_establishDiamonds(establishDiamonds), + m_aaPos(Grid::VertexAttachmentAccessor()), +// m_facDescr(FaceDescriptor()), +// m_volDescr(VolumeDescriptor()), + m_fracInfosBySubset(std::vector()), + //m_sel(Selector()), + m_aAdjMarkerVFP(AttVertFracProp()), + m_aaMarkVrtVFP( Grid::VertexAttachmentAccessor()), +// m_aAdjMarkerVFP(AttVertFracProp()), + m_aaMarkEdgeVFP(Grid::EdgeAttachmentAccessor()), + m_aAdjMarkerB(ABool()), + m_aaMarkFaceB(Grid::FaceAttachmentAccessor()), + m_originalFractureFaces(std::vector()), +// m_attVrtVec(AttVrtVec()), +// m_aaVrtVecVol( Grid::VolumeAttachmentAccessor() ), + m_aAdjInfoEdges(AttVecEdge()), + m_aAdjInfoFaces(AttVecFace()), +// m_aAdjInfoVols(AttVecVol()), + m_aaVrtInfoAssoEdges( Grid::VertexAttachmentAccessor()), + m_aaVrtInfoAssoFaces( Grid::VertexAttachmentAccessor()), +// m_aaVrtInfoAssoVols( Grid::VertexAttachmentAccessor()), + m_aAdjInfoAVVFT( AttVecVertFracTrip() ), + m_aaVrtInfoFraTri(Grid::VertexAttachmentAccessor()), +// m_vrtxFractrQuadrplVec(VrtxFractrQuadrplArte3DVec()) + m_attVrtVec(AttVrtVec()), + m_aaVrtVecVol( Grid::VolumeAttachmentAccessor() ), + m_vecCrossVrtInf(std::vector()), + m_aAdjVolElmInfo(AttVecAttachedVolumeElemInfo()), + m_aaVolElmInfo(Grid::VertexAttachmentAccessor()) +{ +// // Notloesung, nicht in die erste Initialisierung vor geschweifter Klammer, da copy constructor privat + m_sel = Selector(); +} + + +ArteExpandFracs3D::~ArteExpandFracs3D() +{ + // Auto-generated destructor stub +} + +bool ArteExpandFracs3D::run() +{ + if( ! initialize() ) + return false; + + UG_LOG("initialisiert" << std::endl); + + if( ! setSelector() ) + return false; + + UG_LOG("selektiert" << std::endl); + + if( ! attachMarkers() ) + return false; + + UG_LOG("attached" << std::endl); + + if( ! countAndSelectFracBaseNums() ) + return false; + + UG_LOG("gezaehlt" << std::endl); + + + if( ! assignOrigFracInfos() ) + return false; + + UG_LOG("assigniert" << std::endl); + + + if( ! establishNewVrtBase() ) + return false; + + UG_LOG("etabliert" << std::endl); + + if( ! generateVertexInfos() ) + return false; + + UG_LOG("generiert" << std::endl); + + if( ! createConditionForNewVrtcs() ) + return false; + + UG_LOG("kreiert" << std::endl); + + if( ! loop2EstablishNewVertices() ) + return false; + + UG_LOG("loopiert" << std::endl); + + UG_LOG("under construction " << std::endl); + + if( ! createNewElements() ) + return false; + + UG_LOG("new elements created " << std::endl); + + if( ! detachMarkers() ) + return false; + + UG_LOG("detachiert" << std::endl); + + return true; +} + +bool ArteExpandFracs3D::initialize() +{ + UG_LOG("initialize " << std::endl); + + // access position attachment + if(!m_grid.has_vertex_attachment(aPosition) ) + { + UG_LOG("Error in ExpandFractures Arte 3D: Missing position attachment"); + return false; + } + + m_aaPos = Grid::VertexAttachmentAccessor(m_grid, aPosition); + + // make sure that the required options are enabled. + + if( ! m_grid.option_is_enabled(VOLOPT_AUTOGENERATE_FACES) ) + { + UG_LOG("WARNING in Arte 3D init : grid option VOLOPT_AUTOGENERATE_FACES autoenabled.\n"); + m_grid.enable_options(VOLOPT_AUTOGENERATE_FACES); + } + + if( ! m_grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES) ) + { + UG_LOG("WARNING in Arte 3D init: grid option FACEOPT_AUTOGENERATE_EDGES autoenabled.\n"); + m_grid.enable_options(FACEOPT_AUTOGENERATE_EDGES); + } + + // vectors that allow to access fracture properties by subset index + m_fracInfosBySubset = std::vector( m_sh.num_subsets(), FractureInfo(-1, -1, 0) ); + + for( size_t i = 0; i < m_fracInfos.size(); ++i) + { + if( m_fracInfos[i].subsetIndex >= m_sh.num_subsets()) + { + throw(UGError("Bad subsetIndex in given fracInfos.")); + } + + m_fracInfosBySubset[ m_fracInfos[i].subsetIndex] = m_fracInfos[i]; + + } + + + return true; +} + + +bool ArteExpandFracs3D::setSelector() +{ + // Collect surrounding volumes, faces and edges of all fractures in a selector + // and select fracture faces, edges and vertices too. + +// m_sel = Selector(m_grid); + + m_sel.assign_grid(m_grid); + + m_sel.enable_autoselection(false); + m_sel.enable_selection_inheritance(true); //required for select and mark, disabled later + m_sel.enable_strict_inheritance(false); + +// bool strictInherit = m_sel.strict_inheritance_enabled(); +// +// UG_LOG("strikte Inheritenz ist " << strictInherit << " und false ist " << false << std::endl); + + return true; +} + + +bool ArteExpandFracs3D::attachMarkers() +{ + // first part + + // attachment pair boundary is fracture, number fractures crossing + + m_aAdjMarkerVFP = AttVertFracProp(); + +// support::VertexFracturePropertiesVol vfp0; // false, 0 ); + VertxFracPropts vfp0; // false, 0 ); + // default value: no boundary fracture, no fractures crossing + + m_grid.attach_to_vertices_dv( m_aAdjMarkerVFP, vfp0 ); + m_aaMarkVrtVFP = Grid::VertexAttachmentAccessor ( m_grid, m_aAdjMarkerVFP ); + + //m_aAdjMarkerVFP = AttVertFracProp(); + m_grid.attach_to_edges_dv( m_aAdjMarkerVFP, vfp0 ); + + m_aaMarkEdgeVFP = Grid::EdgeAttachmentAccessor( m_grid, m_aAdjMarkerVFP ); + + m_aAdjMarkerB = ABool(); // used to know if an face is frac face + + m_grid.attach_to_faces_dv( m_aAdjMarkerB, false ); + m_aaMarkFaceB = Grid::FaceAttachmentAccessor( m_grid, m_aAdjMarkerB ); + + // second part + +// m_grid.attach_to_volumes(m_attVrtVec); +// m_aaVrtVecVol = Grid::VolumeAttachmentAccessor( m_grid, m_attVrtVec); + + std::vector noEdge; + std::vector noFace; +// std::vector noVol; + + m_aAdjInfoEdges = AttVecEdge(); + m_aAdjInfoFaces = AttVecFace(); +// m_aAdjInfoVols = AttVecVol(); + + m_grid.attach_to_vertices_dv( m_aAdjInfoEdges, noEdge ); + m_aaVrtInfoAssoEdges = Grid::VertexAttachmentAccessor( m_grid, m_aAdjInfoEdges ); + + m_grid.attach_to_vertices_dv( m_aAdjInfoFaces, noFace ); + m_aaVrtInfoAssoFaces = Grid::VertexAttachmentAccessor( m_grid, m_aAdjInfoFaces ); + +// m_grid.attach_to_vertices_dv( m_aAdjInfoVols, noVol ); +// m_aaVrtInfoAssoVols = Grid::VertexAttachmentAccessor( m_grid, m_aAdjInfoVols ); + + + // TODO FIXME + // das fehlt hier , Analogon 2D Fall!!!!!!!!!!!!! der geht hier eigentlich weiter + // die Vertizes, Faces und Edges, die mit einer Kluft zu tun haben + // using VertFracTrip = VertexFractureTriple; + // using VecVertFracTrip = std::vector; + // VecVertFracTrip vertexNoInfo; + + // AttVecVertFracTrip m_aAdjInfoAVVFT; + + VecVertFracTrip vertexNoInfo; + + m_aAdjInfoAVVFT = AttVecVertFracTrip(); + + m_grid.attach_to_vertices_dv( m_aAdjInfoAVVFT, vertexNoInfo ); + + m_aaVrtInfoFraTri = Grid::VertexAttachmentAccessor(m_grid, m_aAdjInfoAVVFT ); + + + // associate a vector of vertices for each volume adjacent to the frac. + // An entry later will contain the new vertex, if the + // corresponding vertex is an inner fracture vertex, and nullptr if not. + + m_attVrtVec = AttVrtVec(); + + m_grid.attach_to_volumes(m_attVrtVec); + + m_aaVrtVecVol = Grid::VolumeAttachmentAccessor(m_grid, m_attVrtVec); + + + VecAttachedVolumeElemInfo noVolInfo; + + m_aAdjVolElmInfo = AttVecAttachedVolumeElemInfo(); + + m_grid.attach_to_vertices_dv(m_aAdjVolElmInfo,noVolInfo); + + m_aaVolElmInfo = Grid::VertexAttachmentAccessor(m_grid, m_aAdjVolElmInfo); + + return true; +} + + +bool ArteExpandFracs3D::detachMarkers() +{ + m_grid.detach_from_vertices( m_aAdjMarkerVFP ); + m_grid.detach_from_edges( m_aAdjMarkerVFP ); + m_grid.detach_from_faces( m_aAdjMarkerB ); + + m_grid.detach_from_vertices( m_aAdjInfoEdges ); + m_grid.detach_from_vertices( m_aAdjInfoFaces ); +// m_grid.detach_from_vertices( m_aAdjInfoVols ); + + m_grid.detach_from_vertices( m_aAdjInfoAVVFT ); + + m_grid.detach_from_volumes( m_attVrtVec ); + + m_grid.detach_from_vertices(m_aAdjVolElmInfo); + + return true; +} + +bool ArteExpandFracs3D::countAndSelectFracBaseNums() +{ + UG_LOG("countandselect" << std::endl); + + for(size_t i_fi = 0; i_fi < m_fracInfos.size(); ++i_fi ) + { + int fracIndSudo = m_fracInfos[i_fi].subsetIndex; + + UG_LOG("sudo ind " << fracIndSudo << std::endl); + + for( FaceIterator iter = m_sh.begin(fracIndSudo); iter != m_sh.end(fracIndSudo); ++iter ) + { + Face* fac = *iter; + +// UG_LOG("Gesicht " << m_aaPos[fac] << std::endl); + + vector3 facCenter = CalculateCenter( fac, m_aaPos ); + UG_LOG("fac center " << facCenter << std::endl); + + for( IndexType i = 0; i < fac->num_vertices(); i++ ) + UG_LOG("Vertex " << i << " -> " << m_aaPos[fac->vertex(i)] << std::endl); + + UG_LOG("alle Vertizes" << std::endl); + + m_sel.select(fac); + + UG_LOG("selektiert msel " << fac << std::endl); + + m_aaMarkFaceB[fac] = true; + + UG_LOG("mark bool " << fac << std::endl); + +// return true; + + std::vector facEdges; + + UG_LOG("kollektiere ecken" << std::endl); + + CollectEdges( facEdges, m_grid, fac ); + + IndexType d_anzahlEcken = facEdges.size(); + + if( d_anzahlEcken == 0 ) + { + UG_LOG("keine Ecken " << std::endl); + return true; + } + + UG_LOG("Anzahl Ecken " << d_anzahlEcken << std::endl); + + for( auto const & edg : facEdges ) + { + m_aaMarkEdgeVFP[edg]++; + + if( IsBoundaryEdge3D( m_grid, edg ) ) + m_aaMarkEdgeVFP[edg].setIsBndFracVertex(); + + m_sel.select(edg); + } + + UG_LOG("Ecken gesammelt " << std::endl); + + for(size_t i = 0; i < fac->num_vertices(); ++i) + { + Vertex* vrt = fac->vertex(i); + m_sel.select(vrt); + // TODO FIXME hier anpassen, herausfinden ob fracture geschlossen + // oder inneres Ende, und Anzahl der umgebenden fractures bestimmen!!! + + auto & vrtxFracPrps = m_aaMarkVrtVFP[ vrt ]; + +// m_aaMarkVrtVFP[vrt]++; + vrtxFracPrps++; + // vielleicht auch in nachfolgendem Loop über alle selektierten Vertizes, + // wo dann die attached faces abgelaufen und dabei die Subdomain Nummer ermittelt wird + // kann eventuell sogar im Hauptloop gemacht werden, dann könnte man praktisch + // das alte vertex fracture properties von 2D weiter verwenden, noch zu klären + +// m_aaMarkVrtVFP[vrt].addFractSudo(fracIndSudo); + vrtxFracPrps.addFractSudo(fracIndSudo); + + if( IsBoundaryVertex3D(m_grid, vrt) ) + { +// m_aaMarkVrtVFP[vrt].setIsBndFracVertex(); + vrtxFracPrps.setIsBndFracVertex(); + } + + // die Ecken heraus filtern, die mit diesem Vertex assoziert sind + std::vector attEdg; + + for( auto const & ae: facEdges ) + { + if( EdgeContains(ae,vrt) ) + attEdg.push_back(ae); + } + + if( attEdg.size() == 2 ) + { + EdgePair edgPr( attEdg[0], attEdg[1] ); + + AttachedFractFaceEdgeSudo afes( fac, edgPr, fracIndSudo ); + + vrtxFracPrps.addAttachedFractElem(afes); + } + else + { + UG_THROW("number of attached edges wrong " << std::endl); + } + + } + } + } + +#if 0 + + // TODO FIXME das ist was komisches, was von Prof. Reiter da ist, es werden edges gesplittet, für was? + + edges.clear(); + for(EdgeIterator iter = sel.begin(); + iter != sel.end(); ++iter) + { + Edge* e = *iter; + if(aaMarkVRT[e->vertex(0)] != 2 && + aaMarkVRT[e->vertex(1)] != 2 && + aaMarkEDGE[e] > 1) + { + edges.push_back(e); + } + } + + for(size_t i = 0; i < edges.size(); ++i){ + vector3 center = CalculateCenter(edges[i], aaPos); + RegularVertex* v = SplitEdge(grid, edges[i], false); + aaPos[v] = center; + aaMarkVRT[v] = 2; + sel.select(v); + // assign adjacency values for associated selected edges (2 to each) + for(Grid::AssociatedEdgeIterator iter = grid.associated_edges_begin(v); + iter != grid.associated_edges_end(v); ++iter) + { + if(sel.is_selected(*iter)) + aaMarkEDGE[*iter] = 2; + } + } + + +#endif + + UG_LOG("neuer Beginn" << std::endl); +// return true; + + // TODO FIXME hier Loop über alle selektierten Vertizes + // darin für jeden Vertex die adjungierten Volumen bestimmen ohne Vorbedingung + // dann den Loop zur Vorbereitung des StammiBene Algorithmus aufrufen + // für jeden Vertex wieder, der hat dann schon die Basisinfo von der + // VecAttachedVolumeElemInfo Klasse an jedem Vertex + // auf die aufbauend fügt er Fracture Manifolds, general manifolds, + // und NEU auch boundary manifolds in einer eigenen Member dazu + // danach wird für jeden Vertex der StammiBene Algorithmus aufgerufen + // später werden Boundary Faces wie eine eigene Subdomain Ebene + // mit Expansion null behandelt + // was an Knicken aussen an boundary zu tun ist, noch zu überlegen + // XXXXXXXXXXXXXXXXXXXXXXXX hier fangen wir an HHHHHHHHHHHHHHHH + + for( VertexIterator iter = m_sel.begin(); iter != m_sel.end(); ++iter) + { + Vertex* vrt = *iter; + + bool wahl = true; + + // TODO FIXME + // hier den Stammi-bene Algorithmus einfügen + // Search the merged manifold interatively between each basic element + // dazu noch für boundary faces ein eigenes member einführen, + // das vom Typ General statt fracture manifold ist, aber sonst wie general + // später wirken die boundary faces wie eine fracture, die aber + // um den Wert null nur verschoben wird + // die Anzahl der Segmente bestimmt, ob der Vertex gewählt wird + // damit er gewählt wird, muss mehr als ein Segment vorhanden sein + + auto & vrtxFracPrps = m_aaMarkVrtVFP[ vrt ]; + +// bool isBnd = m_aaMarkVrtVFP[ vrt ].getIsBndFracVertex(); + bool isBnd = vrtxFracPrps.getIsBndFracVertex(); +// auto numCrosFrac = m_aaMarkVrtVFP[ vrt ].getNumberFracEdgesInVertex(); + + VertxFracPropts::VrtxFracStatus vfpStatus = vrtxFracPrps.getVrtxFracStatus(); + + if( vfpStatus == VertxFracPropts::noFracSuDoAtt ) + UG_THROW("vertex selected and no frac " << std::endl); + + // TODO FIXME das ist richtig für den 2D Fall, aber passt das auch im 3D Fall???? nein, da SudoFrac Zahl nötig +// if( ! isBnd && numCrosFrac == 1 ) +// if( ! isBnd && vfpStatus == VertxFracPropts::oneFracSuDoAtt && ) + // TODO FIXME was, wenn ein Teil geschlossen ist der fractures, und ein anderer nicht??? + //static_assert(false); + + + // bestimmen, ob die vertizes der fracture vertizes von ihrer subdomain komplett umzingelt werden + // muss vor hier geklärt werden!!! + +// VecPairSudoBool sudoIsSourrounded; + + + + bool allClosed = false; + + allClosed = isVrtxSurroundedByFracFaces( vrt, vrtxFracPrps ); //, sudoIsSourrounded ); +// // case boundary: figure out if the vertex is surrounded by frac faces of which two end in +// // boundary edges, similar the case when the boundary face has itself two +// // boundary edges, where the vertex is connected to both of them, then it is easy + +// if( ! allClosed ) +// return false; + +// return true; + +// if( ! isBnd ) +// { +// UG_LOG("test if closed" << std::endl); +// +//// m_sh.assign_subset(vrt,m_sh.num_subsets()); +// +//// UG_LOG("test if closed assign" << std::endl); +// +// +// } +// else +// { +// +// +// +//// allClosed = +// } +// // TODO FIXME auch bei einer boundary muss das gemacht werden! + + UG_LOG("getestet if closed " << m_aaPos[vrt] << std::endl); + + +// return true; + + if( allClosed == vrtxFracPrps.getInfoAllFracturesSameClosedState() ) + UG_THROW("da ist was schief gegangen " << std::endl); + + // TODO FIXME ist das so richtig? kann sein, dass das zu vereinfacht ist!!!!! + // sudo is suourrounded muss übertragen werden TODO FIXME +// if( ! isBnd && vrtxFracPrps.getInfoAllFracturesSameClosedState() ) + // das !isBnd im 2D Fall wichtig, hier bestimmt auch, wie verallgemeinern? + // bei mehreren subdoms eventuell komplizierter, kommt aber hoffentlich am Rand nicht vor...... + if( vrtxFracPrps.getInfoAllFracturesSameClosedState() ) + { + wahl = false; + } + + UG_LOG("SELEKTIERE " << m_aaPos[vrt] << " -> " << vrtxFracPrps.getInfoAllFracturesSameClosedState() << std::endl); + + // was, wenn numCrossFrac == 0 ist? + // wieso werden die boundary vrt ausgeschlossen, oder sollen die nicht ausgeschlossen werden? + // schon im 2D Fall unklar, hier noch wirrer!!! TODO FIXME + + if( wahl ) +// if( m_aaMarkVrtVFP[vrt].getNumberFracEdgesInVertex() > 1 ) // TODO FIXME stimmt das so? + { + // select all associated edges, faces and volumes + m_sel.select( m_grid.associated_edges_begin(vrt), m_grid.associated_edges_end(vrt) ); + m_sel.select( m_grid.associated_faces_begin(vrt), m_grid.associated_faces_end(vrt) ); + m_sel.select( m_grid.associated_volumes_begin(vrt), m_grid.associated_volumes_end(vrt) ); + + std::vector assoEdg; + std::vector assoFac; + std::vector assoVol; + VecAttachedVolumeElemInfo assoVolElemInfo; + + for( std::vector::iterator iterEdg = m_grid.associated_edges_begin(vrt); + iterEdg != m_grid.associated_edges_end(vrt); + iterEdg++ ) + { + assoEdg.push_back(*iterEdg); + } + + for( std::vector::iterator iterFac = m_grid.associated_faces_begin(vrt); + iterFac != m_grid.associated_faces_end(vrt); + iterFac++ ) + { + assoFac.push_back(*iterFac); + } + + // TODO FIXME das nach oben verschieben, wo der Stammi Bene Algo sein soll + // die asso edges und faces brauchen wir vielleicht gar nicht + // bzw asso edges und asso faces können hier bleiben wo gewählt wird + // die assoVolElemInfo wird schon oben erzeugt vor der Wahl + // und dann wird die danach folgende Loop Info + for( std::vector::iterator iterVol = m_grid.associated_volumes_begin(vrt); + iterVol != m_grid.associated_volumes_end(vrt); + iterVol++ ) + { + assoVol.push_back(*iterVol); + + AttachedVolumeElemInfo avei(*iterVol); + + assoVolElemInfo.push_back(avei); + } + + m_aaVrtInfoAssoEdges[vrt] = assoEdg; + m_aaVrtInfoAssoFaces[vrt] = assoFac; +// m_aaVrtInfoAssoVols[vrt] = assoVol; + m_aaVolElmInfo[vrt] = assoVolElemInfo; + + } + } + + UG_LOG("vertex Infos Runde eins fertig " << std::endl); + + + // Voraussetzung FÜR StammiBene Aufrufung + // Stammi-Bene-Vorbereitung + for( VertexIterator iter = m_sel.begin(); iter != m_sel.end(); ++iter) + { + Vertex* vrt = *iter; + +// std::vector & attVol = m_aaVrtInfoAssoVols[vrt]; + + auto & vrtxFracPrps = m_aaMarkVrtVFP[ vrt ]; + + VecAttachedFractFaceEdgeSudo vecAttFacEdgSudo = vrtxFracPrps.getAllAttachedFractElems(); + + auto & vecVolElmInfo = m_aaVolElmInfo[vrt]; + + // TODO FIXME eigentlich eine Dummheit, das auf zu teilen in ein VolElemInfo und ein VrtInfoAssoVols + // denn die InfoAssoVols haben als Info nur die Volumen, die VolElmInfos haben die Volumen + // und noch viel mehr Infos zu den Faces und den Edges.... + // mittelfristig die m_aaVrtInfoAssoVols abschaffen und alles auf die AttachedFullDimElemInfo + // übertragen, dann geht der folgende Loop gleich über den Vektor darüber, bzw. gleichbedeutend + // über m_aaVolElmInfo +// for( auto & vol : attVol ) + for( AttachedVolumeElemInfo & attVolElmInfo : vecVolElmInfo ) + { +// AttachedVolumeElemInfo attVolElmInfo( vol ); + Volume * vol = attVolElmInfo.getFulldimElem(); + + // add those faces which are fracture faces + for( auto & afes : vecAttFacEdgSudo ) + { + attVolElmInfo.addFractManifElem(afes, m_grid); + } + + // add those faces which are NOT fracture faces, assign them arbitraryly subdomain -1 + // to indicate that they are not from the manifold, independent of their subdomain + + // collect all volume faces which incorporate the vertex + + std::vector volFacesContainingVrtx; + + for( IndexType iFac = 0; iFac < vol->num_faces(); iFac++ ) + { + Face * fac = m_grid.get_face(vol,iFac); + + if( FaceContains( fac, vrt ) ) + { + volFacesContainingVrtx.push_back( fac ); + } + } + + for( auto const & fac : volFacesContainingVrtx ) + { + // get the edges of the face connected to the vertex + + std::vector vecEdgesFaceVrtx; + + // need to be two edges always, check + + for( IndexType iEdge = 0; iEdge < fac->num_edges(); iEdge++ ) + { + Edge * edg = m_grid.get_edge(fac,iEdge); + + if( EdgeContains(edg,vrt) ) + { + vecEdgesFaceVrtx.push_back(edg); + } + } + + if( vecEdgesFaceVrtx.size() != 2 ) + { + UG_LOG("edge number Unsinn " << vecEdgesFaceVrtx.size() << std::endl); + UG_THROW("edge number Unsinn " << vecEdgesFaceVrtx.size() << std::endl); + return false; + } + + EdgePair edgesFaceVrtx( vecEdgesFaceVrtx[0], vecEdgesFaceVrtx[1] ); + + // test the subdomain first, if from the subdomains of the cleft manifolds + + IndexType sudoThisFace = m_sh.get_subset_index(fac); + + std::vector const & sudoList = vrtxFracPrps.getSudoList(); + + // test if sudo of face belongs to the fracture face subdom list + + bool belongsToFracFaceSudo = false; + + for( auto const & sudoFrac : sudoList ) + { + if( sudoFrac == sudoThisFace ) + { + belongsToFracFaceSudo = true; + break; + } + } + + + if( belongsToFracFaceSudo ) + { + // if it belongs, construct it again and test if it already belongs to the fracture faces + // MUST be already part of the list, else major error appeared! + + AttachedFractFaceEdgeSudo afesTest( fac, edgesFaceVrtx, sudoThisFace ); + + if( attVolElmInfo.addFractManifElem( afesTest, m_grid ) ) + { + UG_LOG("manifold element already contained!" << std::endl); + UG_THROW("manifold element already contained!" << std::endl); + return false; + } + + // nothing to do, already added before hoffentlich + + } + else + { + // zeitweilig fehlten alle Faces, die keine fractures sind + // die müssen noch irgendwie als nicht-fracture-faces auch dazu gefügt werden + // die sind in den attached volumes schon enthalten, + // Frage: wie prüfen, dass die gar keine fractures sind, die Infos sollten bekannt sein..... + // wichtig auch, dass nur die faces dazu kommen, die den Vertex enthalten!!! + // irgendwas von der Art "nonFractureFaceInfos" oder sowas dazu hängen, mit Info + // ob schon getouched oder noch nicht..... + + + // we construct the attached manifold, given that it is NOT a fracture manifold + + // notwendig auch, dass es eine Markierungsmöglichkeit gibt dafür, ob + // ein face bei der nächsten weiter inneren Runde quasi äussere Begrenzung sein muss + // gilt sowohl für fracture faces, die können das in erster Runde auch sein am Ende + // der Runde, sowie danach nur noch für nicht-fracture-faces + + AttachedGenerFaceEdgeSudo afesAdd( fac, edgesFaceVrtx ); + + attVolElmInfo.addGenerManifElem( afesAdd, m_grid ); + + } + + } + + } + } + + + UG_LOG("vertex Infos Runde eins fertig Volumen auch" << std::endl); + + return true; +} + +// herausfinden für Sudo der frac, ob bezüglich dieser sudo die faces geschlossen sind, oder ob ein Fracture End vorliegt +bool ArteExpandFracs3D::isVrtxSurroundedByFracFaces( Vertex * const & vrt, VertxFracPropts & vrtxFracPrps ) +//, VecPairSudoBool & sudoSurrounded ) +{ + // TODO FIXME wenn an Boundary, dann auch auf closed open unterscheiden - sowohl wenn nur edge an + // boundary, aber auch wenn ein ganzes face dort, noch unklar, was das bedeutet + + // ganz ähnlich wie im 2D Fall, Loopen, im Kreis drehen, kucken, ob wir vom Anfang ans Ende kommen, + // und ob das Ende der edges wieder der Anfang der edges ist, da wir uns auf faces drehen + + // case boundary: figure out if the vertex is surrounded by frac faces of which two end in + // boundary edges, similar the case when the boundary face has itself two + // boundary edges, where the vertex is connected to both of them, then it is easy + + + VecAttachedFractFaceEdgeSudo vafes = vrtxFracPrps.getAllAttachedFractElems(); + + std::vector sudoList = vrtxFracPrps.getSudoList(); + +// for( auto const & sudo : sudoList ) +// { +// std::vector tmpFaces; +// +// CollectFaces( tmpFace, m_grid, vrt ); +// +// } + + // first compare sudo list, if equal + + std::vector sudosTestList; + + std::vector sudoFound( sudoList.size(), false ); + + UG_LOG("sudo list size " << sudoList.size() << std::endl ); + + UG_LOG("vafes list size VA " << vafes.size() << std::endl ); + + + for( auto const & af : vafes ) + { + bool found = false; + + IndexType sudo = af.getSudo(); + + UG_LOG("sudo af " << sudo << std::endl); + + for( IndexType i = 0; i < sudoList.size(); i++ ) + { + UG_LOG("sudo list dusso is " << sudoList[i] << std::endl); + + if( sudo == sudoList[i] ) + { + sudoFound[i] = true; + found = true; + } + } + + + if( ! found ) + UG_THROW("sudo nicht gefunden muss aber da sein " << std::endl); + } + + UG_LOG("alles gefunden " << std::endl); + + + for( auto const & sf: sudoFound ) + { + if( sf == false ) + { + UG_LOG("Falsch" << std::endl); + UG_THROW("sudo not found but must be there " << std::endl); + } + } + + UG_LOG("alles gefunden Test " << std::endl); + + + // sort faces with respect to subdomain - macht das wirklich Sinn, wie umständlich das gemacht wird jetzt? + + // if we arrive here, all sudos found, the entire circle closing can start + + VecPairSudoBool sudoSurrounded; + + for( auto const & sudo : sudoList ) + { + VecAttachedFractFaceEdgeSudo vecAttFacSudo; + + for( auto const & attFac : vafes ) + { + if( sudo == attFac.getSudo() ) + { + UG_LOG("die sudo gefunden " << sudo << std::endl); + vecAttFacSudo.push_back(attFac); + } + } + + VecAttachedFractFaceEdgeSudo vecAttFacSudoSort; + + + + bool isClosed = false; + + bool isBndVrtx = IsBoundaryVertex3D(m_grid,vrt); + + if( ! isBndVrtx ) + { + UG_LOG("No boundary vertex test closed " << m_aaPos[vrt] << std::endl); + + if( vecAttFacSudo.size() == 1 ) + { + // no need for further investigations for inner faces + isClosed = false; + vecAttFacSudoSort = vecAttFacSudo; + } + else + { + isClosed = sortElemCircleIsClosed( vecAttFacSudo, vecAttFacSudoSort ); + } + + } + else // different treatment boundary vertex + { + // figure out start face with a boundary edge, and figure out an eventual additional boundary edge + // if only one face, then check if two attached boundary edges, then true, else false + + if( vecAttFacSudo.size() == 1 ) + { + AttachedFractFaceEdgeSudo & singleEntry = vecAttFacSudo[0]; + + EdgePair faceEdgs = singleEntry.getPairLowElm(); + + if( IsBoundaryEdge3D(m_grid, faceEdgs.first) && IsBoundaryEdge3D(m_grid, faceEdgs.second ) ) + isClosed = true; + + // very simple + vecAttFacSudoSort = vecAttFacSudo; + + } + else + { + // figure out a begin face with a boundary edge and figure out another face with a boundary edge + + Edge * beginEdge = nullptr; + Edge * endEdge = nullptr; + + IndexType startFaceIndx = 0; + IndexType endFaceIndx = 0; + + for( auto const & afs : vecAttFacSudo ) + { + Face * fac = afs.getManifElm(); + + EdgePair edgs = afs.getPairLowElm(); + + Edge * edgOne = edgs.first; + Edge * edgTwo = edgs.second; + + if( beginEdge == nullptr ) + { + if( IsBoundaryEdge3D(m_grid, edgOne) ) + { + beginEdge = edgTwo; + } + else if( IsBoundaryEdge3D(m_grid, edgTwo) ) + { + beginEdge = edgOne; + } + else + { + startFaceIndx++; + } + } + else + { + if( IsBoundaryEdge3D(m_grid, edgOne) ) + { + endEdge = edgOne; + } + else if( IsBoundaryEdge3D(m_grid, edgTwo) ) + { + endEdge = edgTwo; + } + } + + if( endEdge != nullptr ) + break; + + endFaceIndx++; + } + + if( beginEdge == nullptr && endFaceIndx == vecAttFacSudo.size() ) +// || beginEdge == nullptr || endEdge == nullptr ) + { + UG_LOG("keine boundary edges" << std::endl); + + startFaceIndx = -1; + + if( endEdge != nullptr ) + UG_THROW("Ende nicht null, Anfang null " << std::endl); + } + + UG_LOG("Boundary vertex test closed " << m_aaPos[vrt] << std::endl); + +// int d_num = 0; +// for( auto const & afs : vecAttFacSudo ) +// { +// d_num++; +// Face * fac = afs.getManifElm(); +// m_sh.assign_subset(fac, m_sh.num_subsets()); +// } +// UG_LOG("number of surr faces " << d_num << std::endl ); + +// m_sh.assign_subset(beginEdge,m_sh.num_subsets()); +// m_sh.assign_subset(endEdge,m_sh.num_subsets()); +// m_sh.assign_subset(vecAttFacSudo[startFaceIndx].getManifElm(),m_sh.num_subsets()); + + isClosed = sortElemCircleIsClosed( vecAttFacSudo, vecAttFacSudoSort, startFaceIndx, beginEdge, endEdge ); + + } + } + + UG_LOG("-------------------------------" << std::endl); + UG_LOG("is closed " << isClosed << " at " << m_aaPos[vrt] << std::endl); + UG_LOG("-------------------------------" << std::endl); + +// if( isClosed ) +// { +// m_sh.assign_subset(vrt,3); +// } +// else +// { +// m_sh.assign_subset(vrt,4); +// } + + if( vecAttFacSudo.size() != vecAttFacSudoSort.size() ) + { +// return false; + + UG_THROW("Die Sortierung ist komplett schief gegangen " << std::endl); + } + + // DEBUG Zeug, später entfernen!!!!!! +// for( auto const & afss : vecAttFacSudoSort ) +// { +// Face * fac = afss.getManifElm(); +// +// m_sh.assign_subset(fac, m_sh.num_subsets()); +// } + + PairSudoBool ic( sudo, isClosed ); + + sudoSurrounded.push_back( ic ); + + } + + vrtxFracPrps.setInfoAllFractureSudosIfClosed(sudoSurrounded); + + bool allClosed = true; + + for( auto const & ic : sudoSurrounded ) + { + if( ! ic.second ) + allClosed = false; + } + + return allClosed; +} + +bool ArteExpandFracs3D::sortElemCircleIsClosed( VecAttachedFractFaceEdgeSudo const & vecAttFac, + VecAttachedFractFaceEdgeSudo & vecSortedFac, + int startFacIndexUser, +// int endFacIndexUser, +// IndexType startEdgeIndexUser, +// IndexType endEdgeIndexUser +// Face * const & startFacUser, +// Face * const & endFacUser, + Edge * const & startEdgUser, + Edge * const & endEdgUser + ) +{ + + UG_LOG("Schliesstest" << std::endl); + + IndexType originalSize = vecAttFac.size(); + + if( originalSize == 0 ) + { + UG_THROW("zu klein zum sortieren " << std::endl); + return false; + } + else if ( originalSize == 1 ) + { + UG_THROW("should not happen size 1, should have been mentioned before " << std::endl); + } + + UG_LOG("Kopieren zwecks sortieren " << std::endl); + + for( auto const & af : vecAttFac ) + { + UG_LOG("die sudos innen sind " << af.getSudo() << std::endl); + } + + VecAttachedFractFaceEdgeSudo copyVecAttFac = vecAttFac; + + for( auto const & af : copyVecAttFac ) + { + UG_LOG("die sudos kopiert sind " << af.getSudo() << std::endl); + } + + IndexType beginIndx = 0; + + UG_LOG("begin Index " << beginIndx << std::endl); + + bool simpleConnectionTest = false; + + Edge * startEdgeForced = nullptr; + + if( startFacIndexUser >= 0 ) + { + UG_LOG("Veränderung " << startFacIndexUser << std::endl); + beginIndx = startFacIndexUser; + simpleConnectionTest = true; // user hopefully did it + } + else + { + // we need to ensure to start at a fracture which is not in between, in case that circle not closed + // so ensure that all fracture faces have a fracture face at both edges as neighbour, + // in principle this is already sufficient to answer the question which we want to know + // all the rest here is useless playing in principle + + bool broken = false; + + for( IndexType i = 0; i < vecAttFac.size(); i++ ) + { + IndexType firstSideConnected = 0; + IndexType secondSideConnected = 0; + + AttachedFractFaceEdgeSudo afBase = vecAttFac[i]; + + Face * faceBase = afBase.getManifElm(); + EdgePair edgPairBase = afBase.getPairLowElm(); + + Edge * edgeBaseOne = edgPairBase.first; + Edge * edgeBaseTwo = edgPairBase.second; + + for( IndexType j = 0; j < vecAttFac.size(); j++ ) + { + if( i != j ) + { + AttachedFractFaceEdgeSudo afCompr = vecAttFac[j]; + + Face * faceCompr = afCompr.getManifElm(); + EdgePair edgPairCompr = afCompr.getPairLowElm(); + + Edge * edgeComprOne = edgPairCompr.first; + Edge * edgeComprTwo = edgPairCompr.second; + + if( edgeComprOne == edgeBaseOne || edgeComprTwo == edgeBaseOne ) + firstSideConnected++; + + if( edgeComprOne == edgeBaseTwo || edgeComprTwo == edgeBaseTwo ) + secondSideConnected++; + } + + } + + if( vecAttFac.size() > 2 && ( firstSideConnected > 1 || secondSideConnected > 1 ) ) + { + UG_THROW("zu oft verbunden " << std::endl ); + } + else if( vecAttFac.size() == 2 && ( firstSideConnected > 2 || secondSideConnected > 2 ) ) + { + UG_THROW("zu oft verbunden " << std::endl ); + } + else if( firstSideConnected == 0 ) + { + // face is open into one direction, already clear that not closed!! + + beginIndx = i; + simpleConnectionTest = false; + startEdgeForced = edgeBaseTwo; + UG_LOG("forcieren 1 " << std::endl); + broken = true; + break; + } + else if( secondSideConnected == 0 ) + { + // face is open into one direction, already clear that not closed!! + + beginIndx = i; + simpleConnectionTest = false; + startEdgeForced = edgeBaseOne; + UG_LOG("forcieren 2 " << std::endl); + broken = true; + break; + } + else if( firstSideConnected == 1 && secondSideConnected == 1 ) + { + simpleConnectionTest = true; // as long as a look + } + else + { + UG_THROW("komischer Verbindungsfall " << std::endl); + } + + if( broken ) + break; + + } + } + + UG_LOG("begin Index X " << beginIndx << std::endl); + + +// AttachedFractFaceEdgeSudo initialAFES = *(copyAttFac.begin()); + AttachedFractFaceEdgeSudo initialAFES = copyVecAttFac[beginIndx]; + + IndexType sudo = initialAFES.getSudo(); + + UG_LOG("sudo " << beginIndx << " ist " << sudo << std::endl); + + Face * beginFacLoop = initialAFES.getManifElm(); + + // TODO FIXME wird das irgendwo verwendet? wieso nicht? +// Face * endFacLoop = nullptr; + +// if( endFacIndexUser != -1 ) +// { +// endFacLoop = copyVecAttFac[endFacIndexUser].getManifElm(); +// } + + EdgePair beginEdges = initialAFES.getPairLowElm(); + + Edge * beginEdgeLoop = beginEdges.second; + Edge * targetedEndEdgeLoop = beginEdges.first; // should be closed! should end at same edge as it begins! + +// +// return true; + +// if( startEdgeIndexUser != -1 ) +// { +// beginEdgeLoop = +// } + +// if( startFacUser != nullptr ) +// { +// beginFacLoop = startFacUser; +// +// } +// + +// if( ! FaceContains(beginFacLoop,beginEdgeLoop->vertex(0)) ) +// UG_THROW("Sortierung Gesicht hat nicht die gewünschte Ecke " << std::endl); +// +// if( endFacUser != nullptr ) +// { +// endFacLoop = endFacUser; +// } +// + + if( startEdgUser != nullptr ) + { + beginEdgeLoop = startEdgUser; + + // check if part of the begin face! + + if( ! FaceContains( beginFacLoop, beginEdgeLoop ) ) + UG_THROW("Anfangsgesicht hat keine Anfangsecke " << std::endl); + + } + + if( startEdgeForced != nullptr ) + { + beginEdgeLoop = startEdgeForced; + + // check if part of the begin face! + + if( ! FaceContains( beginFacLoop, beginEdgeLoop ) ) + UG_THROW("Anfangsgesicht hat keine Anfangsecke forced " << std::endl); + +// m_sh.assign_subset(startEdgeForced, m_sh.num_subsets()); + UG_LOG("forciert " << std::endl); + +// return false; + + } + + + if( endEdgUser != nullptr ) + { + // check if part of end face at end of loop somehow + targetedEndEdgeLoop = endEdgUser; + + if( startEdgeForced != nullptr ) + UG_THROW("das muss schief gehen, Chaos " << std::endl); + +// if( endFacLoop != nullptr ) +// { +// if( ! FaceContains( endFacLoop, targetedEndEdgeLoop ) ) +// UG_THROW("Endgesicht hat keine Endecke " << std::endl); +// } + } +// + + // DEBUG +// m_sh.assign_subset(beginFacLoop, m_sh.num_subsets()); +// m_sh.assign_subset(beginEdgeLoop, m_sh.num_subsets()); +// m_sh.assign_subset(targetedEndEdgeLoop, m_sh.num_subsets()); + + + // Du musst sortieren. Du musst einen Zeitplan machen. Das kann man lernen. Du kannst das selber machen. + + IndexType countedCrossedFaces = 1; + + vecSortedFac.push_back( initialAFES ); + + copyVecAttFac.erase( copyVecAttFac.begin() + beginIndx ); + +// Face * face2Append = beginFacLoop; + Edge * startEdge2Append = beginEdgeLoop; + +// m_sh.assign_subset(startEdge2Append,m_sh.num_subsets()); + + UG_LOG("while loop anfangen " << std::endl); + + IndexType d_whi = 0; + + while( copyVecAttFac.size() != 0 ) + { + + UG_LOG("in while loop " << d_whi << std::endl); + d_whi++; + + IndexType foundCommEdg = 0; + + Edge * nextEdge = nullptr; + +// for( auto const & caf : copyVecAttFac ) + for( VecAttachedFractFaceEdgeSudo::iterator itAttFES = copyVecAttFac.begin(); + itAttFES != copyVecAttFac.end(); + itAttFES++ + ) + { + AttachedFractFaceEdgeSudo caf = *itAttFES; + + Face * d_Fac = caf.getManifElm(); + +// m_sh.assign_subset(d_Fac,m_sh.num_subsets()); + + EdgePair edgPr = caf.getPairLowElm(); + + Edge * edgOne = edgPr.first; + Edge * edgTwo = edgPr.second; + +// m_sh.assign_subset(edgOne,m_sh.num_subsets()); +// m_sh.assign_subset(edgTwo,m_sh.num_subsets()); + +// return true; + + IndexType hasEdge = 0; + + Edge * overNextEdge = nullptr; + + if( edgOne == startEdge2Append ) + { + nextEdge = edgOne; + overNextEdge = edgTwo; + hasEdge++; + } + + if( edgTwo == startEdge2Append ) + { + nextEdge = edgTwo; + overNextEdge = edgOne; + hasEdge++; + } + + if( hasEdge > 1 ) + UG_THROW("zu viele Ecken und Kanten " << std::endl); + + if( hasEdge == 1 ) + { + Face * fac2App = caf.getManifElm(); +// m_sh.assign_subset(fac2App,m_sh.num_subsets()); + EdgePair edgesNextFace( edgOne, edgTwo ); + AttachedFractFaceEdgeSudo nextAttFES( fac2App, edgesNextFace, sudo ); + + vecSortedFac.push_back(nextAttFES); + + copyVecAttFac.erase(itAttFES); + foundCommEdg++; + startEdge2Append = overNextEdge; + + break; + } + + + } + + if( foundCommEdg > 1 ) + UG_THROW("Kein Anschluss unter dieser Nummer " << std::endl); + + if( foundCommEdg == 0 ) + { + UG_LOG("Kreislauf nicht geschlossen " << std::endl); + + if( nextEdge != nullptr ) + UG_THROW("nicht konsistent, null und null " << std::endl); + + break; + } + + if( nextEdge == nullptr ) + { + if( foundCommEdg != 0 ) + UG_THROW("nicht konsistent, null und null v2 " << foundCommEdg << " -> " << nextEdge << std::endl); + + // return false; + } + + } + + if( originalSize != vecSortedFac.size() ) + { + UG_THROW("Sortierung hat nicht funktioniert " << std::endl); + return false; + } + + if( startEdge2Append != targetedEndEdgeLoop ) + { + if( simpleConnectionTest ) + UG_THROW("obwohl offen oder vorgegeben trotzdem Ziel nicht erreicht?" << std::endl); + + UG_LOG("Ende nicht erreicht, Kreis nicht geschlossen, innerer Rand vermutlich" << std::endl); + return false; + } + + + UG_LOG("Kreislauf Faces 3D Test ob Knoten umrandet geschlossen " << std::endl); + + return true; +} + +bool ArteExpandFracs3D::assignOrigFracInfos() +{ + m_originalFractureFaces.clear(); + + for( FaceIterator iter = m_sel.begin(); iter != m_sel.end(); ++iter) + { + if( m_aaMarkFaceB[*iter] == true ) + m_originalFractureFaces.push_back(*iter); + } + + m_fracInfosBySubset = std::vector( m_sh.num_subsets(), FractureInfo(-1, -1, 0) ); + + for(size_t i = 0; i < m_fracInfos.size(); ++i) + { + if( m_fracInfos[i].subsetIndex >= m_sh.num_subsets()) + { + throw(UGError("Bad subsetIndex in given fracInfos.")); + } + + m_fracInfosBySubset[ m_fracInfos[i].subsetIndex ] = m_fracInfos[i]; + } + +// disable selection inheritance to avoid infinite recursion. + m_sel.enable_selection_inheritance(false); + + return true; +} + +bool ArteExpandFracs3D::establishNewVrtBase() +{ + // iterate over all surrounding volumes to enable shifted vertices, this loop taken from SR but shortened + + for( VolumeIterator iterSurrVol = m_sel.volumes_begin(); iterSurrVol != m_sel.volumes_end(); ++ iterSurrVol ) + { + Volume* sv = *iterSurrVol; + + std::vector & newVrts = m_aaVrtVecVol[sv]; + newVrts.resize(sv->num_vertices()); + + for(size_t iVrt = 0; iVrt < sv->num_vertices(); ++ iVrt ) + { + newVrts[iVrt] = nullptr; + } + + // erstmal so tun, als ob keine neuen Vertizes erzeugt werden an den alten Vertizes + + } + + return true; +} + +// Analogon zu VertrexFractureInfo in 2D, wo jeder Vertex eine Liste bekommt, wo alle die ihm angehängten +// Ecken, Faces und Volumen gespeichert werden; dazu die Normalen, und vielleicht noch weitere Infos +bool ArteExpandFracs3D::generateVertexInfos() +{ + // TODO FIXME das wird benötigt + + // sowas von der Art als attachement bei den attachments, und dann mit Leben füllen für jeden Vertex + // in dieser Funktion; + // vielleicht braucht es auch Edge Infos, oder nur Edge Infos? +// VecVertFracTrip vertexNoInfo; +// using AttVecVertFracTrip = Attachment; +// AttVecVertFracTrip aAdjInfoAVVFT; +// grid.attach_to_vertices_dv( aAdjInfoAVVFT, vertexNoInfo ); +// Grid::VertexAttachmentAccessor aaVrtInfoFraTri(grid, aAdjInfoAVVFT ); + + // Lebendigmachung in: + // for( auto & fsfpmv : fracSubdom_facePerpendMinVal ) ..... + // von dort lernen!!!!! + + // notwendig: face, normal, volume, edge + + // TODO FIXME das wollen wir nicht, sondern das alte Vertex Fracture Triple +// m_vrtxFractrQuadrplVec = VrtxFractrQuadrplArte3DVec(); + // TODO FIXME diese Einträge erzeugen + + // TODO FIXME kann vielleicht vereinfacht werden durch einen Loop über alle Vertizes, + // und das Abfragen der dort schon gespeicherten vertex property Geschichten, sonst + // wird manches doppelt gemoppelt + + for( auto const & fracInf : m_fracInfos ) + { + IndexType fracSudo = fracInf.subsetIndex; + + +// for(EdgeIterator iterEdg = m_sh.begin(fracIndSudo); iterEdg != m_sh.end(fracIndSudo); iterEdg++ ) +// + for( FaceIterator iterFac = m_sh.begin(fracSudo); iterFac != m_sh.end(fracSudo); iterFac++ ) + { + +// VrtxFractrQuadrplArte3D vrtxFractrQuadrpl; + + Face* fac = *iterFac; + + auto sudoFacInnerLoop = m_sh.get_subset_index(fac); + + if( sudoFacInnerLoop != fracSudo ) + UG_THROW("Subdomain Index Fehler 3D " << std::endl); + +// std::vector verticesFac; +// +// for( IndexType i = 0; i < fac->num_vertices(); i++ ) +// { +// verticesFac.push_back( fac->vertex(i) ); +// } + + + std::vector assoVols; + + // wo kam denn der Käse her? +// if( ! m_grid.option_is_enabled(FACEOPT_STORE_ASSOCIATED_VOLUMES) ) +// UG_THROW("How to collect asso vols?" << std::endl); + + + if(! m_grid.option_is_enabled(VOLOPT_AUTOGENERATE_FACES) ) + { + UG_LOG("WARNING grid option VOLOPT_AUTOGENERATE_FACES autoenabled.\n"); + m_grid.enable_options(VOLOPT_AUTOGENERATE_FACES); + } + + +// for( Grid::AssociatedVolumeIterator volIt = m_grid.associated_volumes_begin(fac); +// volIt != m_grid.associated_volumes_end(fac); +// volIt++ +// ) +// { +// assoVols.push_back(*volIt); +// } + + CollectVolumes(assoVols, m_grid, fac ); +// +// for( auto const & av : assoVols ) +// { +// m_sh.assign_subset(av, m_sh.num_subsets()); +// } +// +// return true; + +// using VolNormPair = std::pair< Volume*, vector3 >; +// +// using VecVolNormPair = std::vector; +// +// VecVolNormPair vecNormalsAwayVol; + +// UG_LOG("Center Face " << CalculateCenter(fac,m_aaPos) << std::endl); + + for( auto const & kuhVol : assoVols ) + { + bool facFound = false; + IndexType numFd = 0; + + for( IndexType iSide = 0; iSide < kuhVol->num_sides(); iSide++ ) + { + Face * kuhFac = m_grid.get_side(kuhVol, iSide); + +// UG_LOG("Center Kuh Face " << CalculateCenter(kuhFac,m_aaPos) << std::endl); + + // TODO FIXME eigentliches Ziel ist es, den face descriptor des Volumens zu finden, + // das mit dem face übereinstimmt, alle anderen Seiten des Volumens sind egal + // Funktion suchen, die ausgehend von einem Face, das ein Volumen begrenzt, + // den zum Volumen passenden FaceDescriptor findet, also auch richtige Orientierung + // der Vertices beinhaltet + // FRAGE TODO FIXME ist ein face descriptor von der Orientierung zum Volumen abhängig + // oder hängt der nur vom Face ab, das eine vorgegebene Oriertierung hat? + + bool checkCoincide = checkIfFacesVerticesCoincide( kuhFac, fac ); + + if( checkCoincide ) + { + facFound = true; + numFd++; + + if( kuhFac != fac ) + UG_LOG("Kuh Fac ist nicht fac " << std::endl); + + FaceDescriptor facDescr; + + // TODO FIXME testen, ob der Face Descriptor von der Orientierung abhängt + // also testen, ob sich der face descriptor ändert, wenn das Volumen + // auf der einen und auf der anderen Seite des faces hängt + // deswegen auch die ganze Prozedur mit den kuhFacs, die hoffentlich + // je nach Volumen anders orientiert sind als das eigentliche Face, + // aber dieselben Vertices haben, also geometrisch gleich sind, aber anders orientiert!!!! + + // TODO FIXME andere Hergehensweise vielleicht: + // von m_aaVrtInfoAssoVols ausgehen, darüber loopen, oder die in einen Vektor stecken, + // wo die Vertices dabei sind, dann kann man sich vielelicht ein paar Klimmzüge sparen, + // vielleicht aber auch nicht..... + + kuhVol->face_desc( iSide, facDescr ); + + vector3 normal; + + CalculateNormal( normal, & facDescr, m_aaPos ); + + vector3 facCenter = CalculateCenter( kuhFac, m_aaPos ); + vector3 kuhCenter = CalculateCenter( fac, m_aaPos ); + vector3 kuhVolCenter = CalculateCenter( kuhVol, m_aaPos); + +// UG_LOG("Normale zum face descriptor " << normal << " , " << facCenter << std::endl); +// UG_LOG("Normale zum Kuhh descriptor " << normal << " , " << kuhCenter << std::endl); +// UG_LOG("Zentrum des Vol") + +// UG_LOG("fac " << fac << std::endl ); +// UG_LOG("kuh " << kuhFac << std::endl ); + + UG_LOG("Normale ist " << normal << " fac " << facCenter + << " vol " << kuhVolCenter << std::endl); + + +// VolNormPair normalsAwayVol( kuhVol, normal ); +// +// vecNormalsAwayVol.push_back( normalsAwayVol ); + + std::vector facEdgs; + + CollectEdges( facEdgs, m_grid, fac); + + for( IndexType iF = 0; iF < fac->num_vertices(); iF++ ) + { + Vertex * vrt = fac->vertex(iF); + + // verticesFac.push_back(vrt); + + std::vector edgOfVrtx; + + for( auto const & edg : facEdgs ) + { + if( EdgeContains(edg, vrt) ) + { + edgOfVrtx.push_back(edg); + } + } + + if( edgOfVrtx.size() == 2 ) + { + EdgePair commonEdges(edgOfVrtx[0], edgOfVrtx[1]); // fill values + // edges commun between face and volume, with the vertex included as well, i.e. two possibilities + + VertFracTrip infoVerticesThisFace( fac, fracSudo, kuhVol, normal, commonEdges ); + + // TODO FIXME hier irgendwie graphische Ausgabe von irgendwas + + m_aaVrtInfoFraTri[vrt].push_back( infoVerticesThisFace ); + + // DEBUG OUTPUT; REMOVE LATER +// m_sh.assign_subset(kuhVol,m_sh.num_subsets()); + } + else + { + UG_THROW("Mein Face das hat keine Ecken, keine Ecken hat mein Face" << std::endl); + } + } + + } + } + + if( ! facFound || numFd != 1 ) + { + UG_THROW("Kein Kuh Volumen gefunden" << std::endl); + return false; + } + + } + + } + + + } + + return true; +} + +bool ArteExpandFracs3D::checkIfFacesVerticesCoincide( Face * const & facOne, Face * const & facTwo ) +{ + + if( facOne->size() != facTwo->size() ) + return false; + + std::vector facOneVrtcs, facTwoVrtcs; + + collectFaceVertices( facOneVrtcs, facOne ); + collectFaceVertices( facTwoVrtcs, facTwo ); + + for( auto const & vrtOne : facOneVrtcs ) + { + bool found = false; + + IndexType numFd = 0; + + for( auto const & vrtTwo : facTwoVrtcs ) + { + if( vrtOne == vrtTwo ) + { + found = true; + numFd++; + } + } + + if( ! found || numFd != 1 ) + return false; + } + + return true; +} + +bool ArteExpandFracs3D::collectFaceVertices( std::vector & facVrt, Face * const & fac ) +{ + if( fac == nullptr ) + return false; + + facVrt.clear(); + + for( IndexType iF = 0; iF < fac->num_vertices(); iF++ ) + { + Vertex * vrt = fac->vertex(iF); + + facVrt.push_back( vrt ); + } + + return true; +} + + + + +// major function of new grid generation, in Keil Style, but functional grid, only the diamonds have to be +// established in additional functionalities independent of this function +bool ArteExpandFracs3D::loop2EstablishNewVertices() +{ + m_vecCrossVrtInf = std::vector(); + + + // TODO FIXME sowas von der Art wird nötig sein als Vorbereitung für die Diamanten, + // Infos darin speichern, vielleicht auch noch notwendig, die Kanten zu speichern oder die faces, + // zu klären im Laufe der Implementation + // std::vector vecCrossVrtInf; + + // zentraler Loop + + // TODO FIXME vorher noch den attVrtVec und den aaVrtVecFac analog implementieren, das fehlt noch!!! + // ebenso seine Befüllung, braucht noch eine Funktion dazwischen, die attachments selber in die + // attach Funktion natürlich + + int untilVrt = 0; + + for( VertexIterator iterV = m_sel.begin(); iterV != m_sel.end(); ++ iterV ) + { + Vertex * oldVrt = *iterV; + + // Position dieses Vertex + vector3 posOldVrt = m_aaPos[oldVrt]; + + // TODO FIXME diese Funktion mit Leben und Analytischer Geometrie 13. Klasse füllen + +// VecVertFracTrip & vecVertFracTrip = m_aaVrtInfoFraTri[oldVrt]; +// +// std::vector & allAssoEdges = m_aaVrtInfoAssoEdges[oldVrt]; +// std::vector & allAssoFaces = m_aaVrtInfoAssoFaces[oldVrt]; +// std::vector & allAssoVolumes = m_aaVrtInfoAssoVols[oldVrt]; + + UG_LOG("vertex at " << posOldVrt << std::endl); + + auto & vrtxFracPrps = m_aaMarkVrtVFP[ oldVrt ]; + + bool vrtxIsBndVrt = vrtxFracPrps.getIsBndFracVertex(); + + UG_LOG("is bndry " << vrtxIsBndVrt << std::endl); + + VertxFracPropts::VrtxFracStatus statusThisVrtx = vrtxFracPrps.getVrtxFracStatus(); + + if( ! vrtxIsBndVrt ) + { + if( vrtxFracPrps.getInfoAllFracturesSameClosedState() ) + { + // gar nix tun, alle offen, innerer Vertex, darf man hier ankommen? NEIN TODO FIXME +// UG_THROW("hier sollten wir nicht angekommen sein " << std::endl); + UG_LOG("vertex nicht geschlossen alle subdoms, nix tun " << m_aaPos[oldVrt] << std::endl); + } + // TODO FIXME: was, wenn ein Zwischending, Mischung? + else if( statusThisVrtx == VrtxFracProptsStatus::noFracSuDoAtt ) + { + UG_THROW("gar keine Frac darf hier nicht ankommen " << std::endl ); + } + else if( statusThisVrtx == VrtxFracProptsStatus::oneFracSuDoAtt ) + { + // TODO FIXME erster Fall, eine Fracture, innen, geschlossen, kann eigentlich nur hier ankommen + UG_LOG("aktuelles Ziel eine sudo ausdehen " << m_aaPos[oldVrt] << std::endl); + + constexpr bool applyGeneralSegmentOrdering = true; + + establishNewVertices( oldVrt ); + +// if( untilVrt == 10 ) +// return false; +// +// untilVrt++; + +// return false; + +// if( applyGeneralSegmentOrdering ) +// { +// UG_LOG("restrict to Tetrahedra, under construction" << std::endl); +// // Zustand Maulbronn, teilweise funktionierend für Hexahedra, aber nicht zuverlässig +// // und nicht funktional für Tetrahedra +// +// } +// else +// { +// UG_LOG("restrict to Hexahedra, works only in part" << std::endl); +// // Zustand Maulbronn, teilweise funktionierend für Hexahedra, aber nicht zuverlässig +// // und nicht funktional für Tetrahedra +// establishNewVertices( oldVrt ); +// // sufficient to tell the vertex, as the attachements are class members and can be asked in the function +// //m_sh.assign_subset(oldVrt,m_sh.num_subsets()); +// } +// else +// { +// UG_THROW("keine Ahnung auf was beschraenkt" << std::endl); +// } + } + else if( statusThisVrtx == VrtxFracProptsStatus::twoFracSuDoAtt ) + { + UG_LOG("zwei sudos kreuzen in " << m_aaPos[oldVrt] << std::endl); + } + else if( statusThisVrtx == VrtxFracProptsStatus::threeFracSuDoAtt ) + { + UG_LOG("drei sudos kreuzen in " << m_aaPos[oldVrt] << std::endl); + + } + else + { + UG_THROW("was für ein Knoten Status???? " << std::endl); + } + } + else // boundary vertex + { + UG_LOG("Boundary vertizes ausdehen folgend erst " << m_aaPos[oldVrt] << std::endl); + } + + // TODO FIXME wichtig: surrounded status closed / open TODO FIXME, sowie Anzahl der schneidenden Fracs + + } + + return false; + + return true; +} + +//////////////////////////////////////////////////////////////////// + +//template <> +//bool ArteExpandFracs3D::establishNewVertices< Tetrahedron, +// ArteExpandFracs3D::VrtxFracProptsStatus::oneFracSuDoAtt +// >( Vertex * const & oldVrt ) +//template< bool APPLY_GENERAL_SEGMENT_ORDERING, +// ArteExpandFracs3D::VrtxFracProptsStatus vfp, +// typename std::enable_if< std::integral_constant> +// > +template <> +bool ArteExpandFracs3D::establishNewVertices< true, + ArteExpandFracs3D::VrtxFracProptsStatus::oneFracSuDoAtt + >( Vertex * const & oldVrt ) +{ + UG_LOG("under construction Tetrahedra limited" << std::endl); + + VecVertFracTrip const & vecVertFracTrip = m_aaVrtInfoFraTri[oldVrt]; + + VecAttachedVolumeElemInfo const & vecAttVolElemInfo = m_aaVolElmInfo[oldVrt]; + + VecAttachedVolumeElemInfo vecAttVolElemInfoCop = vecAttVolElemInfo; // echte KOPIE + + VecAttachedVolumeElemInfo reconstructedVecAttVolElmInf; + + VecSegmentVolElmInfo vecSegVolElmInfo; + + /* + * While Schleifen aufbauen für den + * Search the adjacent surface interatively - Algorithmus + * (Stasi Algorithmus) + * + */ + + IndexType d_segmenteErledigt = 0; + + while( vecAttVolElemInfoCop.size() != 0 ) + { + SegmentVolElmInfo segmentAVEI; + + AttachedVolumeElemInfo & startVolInfoThisSegment = vecAttVolElemInfoCop[0]; + + startVolInfoThisSegment.markIt(); + + Volume * volSta = startVolInfoThisSegment.getFulldimElem(); + + vector3 center; + + if( volSta != nullptr ) + center = CalculateCenter(volSta,m_aaPos); + +// UG_LOG("volume center " << center << std::endl ); + + int d_loopsDone = 0; + + while( vecAttVolElemInfoCop.size() != 0 ) + { + // count number of marked elements + IndexType numMarkedElems = 0; + IndexType markPoint = 0; + +// IndexType lastMarkPt = 0; + IndexType startIndexInner = 0; + + for( AttachedVolumeElemInfo const & volElInfCop : vecAttVolElemInfoCop ) + { + if( volElInfCop.isMarked() ) + { + Volume * vol = volElInfCop.getFulldimElem(); +// m_sh.assign_subset(vol, m_sh.num_subsets()); + + vector3 center = CalculateCenter(vol,m_aaPos); + +// UG_LOG("DAS ZENTRUM " << numMarkedElems << " -> " << center << std::endl); + + startIndexInner = markPoint; + numMarkedElems++; + + } + + markPoint++; + } + + UG_LOG("LOOPS DONE " << numMarkedElems << std::endl); + + if( numMarkedElems == 0 ) + break; + +// int startIndexInner = -1; +// +// for( int i = 0; i < vecAttVolElemInfoCop.size(); i++ ) +// { +// AttachedVolumeElemInfo vi = vecAttVolElemInfoCop[i]; +// +// Volume * vo = vi.getFulldimElem(); +// +// vector3 center = CalculateCenter(vo,m_aaPos); +// +// UG_LOG("DAS ZENTRUM ZAHL VOR " << i << " -> " << center << std::endl); +// +// if( vi.isMarked() ) +// startIndexInner = i; +// } +// +// if( startIndexInner < 0 ) +// { +// UG_THROW("kein Anfang gefunden " << std::endl); +// } +// +//#if 0 +// IndexType startIndexInner = markPoint - 1; +//#endif + AttachedVolumeElemInfo startVolInfoMarkLoop = vecAttVolElemInfoCop[startIndexInner]; + + Volume * stattVoll = startVolInfoMarkLoop.getFulldimElem(); + + vector3 centerX = CalculateCenter(stattVoll,m_aaPos); + + UG_LOG("DAS ZENTRUM DANACH " << startIndexInner << " -> " << centerX << std::endl); + +// m_sh.assign_subset(stattVoll, m_sh.num_subsets()); +#if 0 + for( int i = 0; i < vecAttVolElemInfoCop.size(); i++ ) + { + AttachedVolumeElemInfo vi = vecAttVolElemInfoCop[i]; + + Volume * vo = vi.getFulldimElem(); + + vector3 center = CalculateCenter(vo,m_aaPos); + + UG_LOG("DAS ZENTRUM ZAHL " << i << " -> " << center << std::endl); + + } +#endif + for( AttachedVolumeElemInfo const & possibleOrigVolInfo : vecAttVolElemInfo ) + { + if( possibleOrigVolInfo.hasSameFulldimElem( startVolInfoMarkLoop ) ) + { + segmentAVEI.push_back(possibleOrigVolInfo); + reconstructedVecAttVolElmInf.push_back(possibleOrigVolInfo); + break; + } + } + + vecAttVolElemInfoCop.erase( vecAttVolElemInfoCop.begin() + startIndexInner ); + +// if( d_loopsDone == 1 ) +// return false; + + for( VecAttachedVolumeElemInfo::iterator aveiIt = vecAttVolElemInfoCop.begin(); + aveiIt < vecAttVolElemInfoCop.end(); + aveiIt++ + ) + { + AttachedVolumeElemInfo & possibleNeighbour = *aveiIt; + + if( possibleNeighbour.hasSameFulldimElem( startVolInfoMarkLoop ) ) + { + continue; + } + else + { + bool neighbourFound = possibleNeighbour.testFullDimElmNeighbour( startVolInfoMarkLoop ); + + if( neighbourFound ) + { + Volume * vol = possibleNeighbour.getFulldimElem(); + +// m_sh.assign_subset(vol, m_sh.num_subsets()); + + } + } + } + + + d_loopsDone++; + + + } + + vecSegVolElmInfo.push_back(segmentAVEI); + +// d_segmenteErledigt++; +// +// if( d_segmenteErledigt == 1 ) +// return false; + } + + if( reconstructedVecAttVolElmInf.size() != vecAttVolElemInfo.size() ) + { + UG_LOG("Rekonstruktion schief gegangen " << std::endl); + UG_THROW("Rekonstruktion schief gegangen " << std::endl); + return false; + } + + for( SegmentVolElmInfo const & svei : vecSegVolElmInfo ) + { + // TODO FIXME das hier wieder entfernen, die Subdomain Zuweisung, nur für debug Zwecke + IndexType sudoMax = m_sh.num_subsets(); + + for( AttachedVolumeElemInfo const & vei : svei ) + { + Volume * vol = vei.getFulldimElem(); + + m_sh.assign_subset( vol, sudoMax ); + } + } + + return true; +} + + +//////////////////////////////////////////////////////////////////// + + + +//template <> +//bool ArteExpandFracs3D::establishNewVertices< Hexahedron, +// ArteExpandFracs3D::VrtxFracProptsStatus::oneFracSuDoAtt +// >( Vertex * const & oldVrt ) +//template< bool APPLY_GENERAL_SEGMENT_ORDERING, +// ArteExpandFracs3D::VrtxFracProptsStatus vfp, +// typename std::enable_if< std::integral_constant> +// > +template <> +bool ArteExpandFracs3D::establishNewVertices< false, + ArteExpandFracs3D::VrtxFracProptsStatus::oneFracSuDoAtt + >( Vertex * const & oldVrt ) +{ + VecVertFracTrip & vecVertFracTrip = m_aaVrtInfoFraTri[oldVrt]; + +// std::vector & allAssoEdges = m_aaVrtInfoAssoEdges[oldVrt]; +// std::vector & allAssoFaces = m_aaVrtInfoAssoFaces[oldVrt]; +// std::vector & allAssoVolumes = m_aaVrtInfoAssoVols[oldVrt]; + + // TODO FIXME works if at all only for very simple geometries + // and works only in particular if all asso volumes are part of the triple, i.e. have all a common face with the fracture + // this is not selbstverständlich at all!!!! + + VecVertFracTrip vecVertFracTripCopy = m_aaVrtInfoFraTri[oldVrt]; // copy, not reference! + + + VecVertFracTrip firstSegment; + VecVertFracTrip secondSegment; + + IndexType beginIndex = 0; + + VertFracTrip startTrip = vecVertFracTripCopy[beginIndex]; + + firstSegment.push_back( startTrip ); + + vector3 normalOne = startTrip.getNormal(); + vector3 normalTwo; + // vermute einfach umgekehrte Normale, ein Trick, der im einfachsten Fall geht...... + VecScale(normalTwo,normalOne,-1); + + vecVertFracTripCopy.erase( vecVertFracTripCopy.begin() + beginIndex ); + +// for( VecVertFracTrip::iterator itVFT = vecVertFracTripCopy.begin(); +// itVFT != vecVertFracTripCopy.end(); +// itVFT++ +// ) +// { +// VertFracTrip actualVFT = *itVFT; +// +// vector3 volCenter = CalculateCenter(actualVFT.getFullElm(),m_aaPos); +// UG_LOG("Vol center" << volCenter << std::endl); +// vecVertFracTripCopy.erase(itVFT); +// } +// +// return true; + +// while( vecVertFracTripCopy.size() != 0 ) +// { +// for( VecVertFracTrip::iterator itVFT = vecVertFracTripCopy.begin(); +// itVFT != vecVertFracTripCopy.end(); +// itVFT++ +// ) + for( auto const & actualVFT : vecVertFracTripCopy ) + { +// VertFracTrip actualVFT = *itVFT; + + vector3 normalActual = actualVFT.getNormal(); + + vector3 volCenter = CalculateCenter(actualVFT.getFullElm(),m_aaPos); + UG_LOG("Vol center" << volCenter << std::endl); + + // test direction + + number cosinus2One = VecDot(normalOne,normalActual); + number cosinus2Two = VecDot(normalTwo,normalActual); + + UG_LOG("normal eins " << normalOne << std::endl); + UG_LOG("normal zwei " << normalTwo << std::endl); + UG_LOG("normal actu " << normalActual << std::endl); + + UG_LOG("cosi one " << cosinus2One << std::endl ); + UG_LOG("cosi two " << cosinus2Two << std::endl ); + + + // if cosinus > 0, assume same side + + if( ( cosinus2One >= 0 && cosinus2Two >= 0 ) || ( cosinus2One <= 0 && cosinus2Two <= 0 ) ) + UG_THROW("kann nicht auf zwei Seiten hinken" << std::endl); + + if( cosinus2One >= 0 ) + { + firstSegment.push_back( actualVFT ); + } + else if( cosinus2Two >= 0 ) + { + secondSegment.push_back( actualVFT ); + } + else + { + UG_THROW("muss wo dazu gehoeren wohl" << std::endl); + } + } +// vecVertFracTripCopy.erase(itVFT); +// } +// } + + // computer averaged normal + + vector3 normalsOneSummed(0,0,0); + vector3 normalsTwoSummed(0,0,0); + + + for( auto const & seg: firstSegment ) + { + vector3 tmpVec = normalsOneSummed; + VecAdd(normalsOneSummed,tmpVec,seg.getNormal()); + } + + for( auto const & seg: secondSegment ) + { + vector3 tmpVec = normalsTwoSummed; + VecAdd(normalsTwoSummed,tmpVec,seg.getNormal()); + } + + vector3 normalsOneAveraged; + vector3 normalsTwoAveraged; + + if( firstSegment.size() != 0 ) + { + VecScale(normalsOneAveraged, normalsOneSummed, 1./firstSegment.size()); + } + + if( secondSegment.size() != 0 ) + { + VecScale(normalsTwoAveraged, normalsTwoSummed, 1./secondSegment.size()); + } + + // get to know width of fracture + + IndexType suse = startTrip.getSudoElm(); + + number width = m_fracInfosBySubset[suse].width; + + number scal = width / 2.; + + vector3 scaledNormalOne, scaledNormalTwo; + VecScale( scaledNormalOne, normalOne, - scal ); + VecScale( scaledNormalTwo, normalTwo, - scal ); + // Minuszeichen wichtig, sonst wird in die falsche Richtung gedrückt, und die Volumen gehen über die fracture + // raus und werden grösser, anstatt kleiner zu werden..... + + vector3 posOldVrt = m_aaPos[oldVrt]; + + vector3 posNewVrtOne, posNewVrtTwo; + + VecAdd( posNewVrtOne, posOldVrt, scaledNormalOne); + VecAdd( posNewVrtTwo, posOldVrt, scaledNormalTwo); + + Vertex * newShiftVrtxOne = *m_grid.create(); + Vertex * newShiftVrtxTwo = *m_grid.create(); + + m_aaPos[newShiftVrtxOne] = posNewVrtOne; + m_aaPos[newShiftVrtxTwo] = posNewVrtTwo; + + UG_LOG("Created new vertex 1 at " <& newVrts4Fac = m_aaVrtVecVol[ vol ]; + + for(size_t indVrt = 0; indVrt < (vol)->num_vertices(); indVrt++ ) + { + Vertex* volVrt = (vol)->vertex(indVrt); + + if( volVrt == oldVrt ) + { + newVrts4Fac[ indVrt ] = newShiftVrtxOne; + } + } + } + + for( auto const & ses : secondSegment ) + { + Volume * vol = ses.getFullElm(); + + std::vector& newVrts4Fac = m_aaVrtVecVol[ vol ]; + + for(size_t indVrt = 0; indVrt < (vol)->num_vertices(); indVrt++ ) + { + Vertex* volVrt = (vol)->vertex(indVrt); + + if( volVrt == oldVrt ) + { + newVrts4Fac[ indVrt ] = newShiftVrtxTwo; + } + } + } + + + return true; +} + +//////////////////////////////////////////////////////////////////// + + +bool ArteExpandFracs3D::createConditionForNewVrtcs() +{ + + // iterate over all surrounding volumes to enable volume changes, this loop taken from SR but shortened + for(VolumeIterator iterSurrVol = m_sel.volumes_begin(); iterSurrVol != m_sel.volumes_end(); iterSurrVol++ ) + { + Volume * sv = *iterSurrVol; + + std::vector& newVrts = m_aaVrtVecVol[sv]; + newVrts.resize(sv->num_vertices()); + + for(size_t iVrt = 0; iVrt < sv->num_vertices(); iVrt++ ) + { + newVrts[iVrt] = nullptr; + } + // erstmal so tun, als ob keine neuen Vertizes erzeugt werden an den alten Vertizes + } + + + return true; +} + +///////////////////////////////////////////////////////////// + +bool ArteExpandFracs3D::createNewElements() +{ + // practically copied from Sebastian, as this concept is fine + + // create new elements + + // holds local side vertex indices + std::vector locVrtInds; + + // first we create new edges from selected ones which are connected to + // inner vertices. This allows to preserve old subsets. + // Since we have to make sure that we use the right vertices, + // we have to iterate over the selected volumes and perform all actions on the edges + // of those volumes. + for(VolumeIterator iter_sv = m_sel.volumes_begin(); iter_sv != m_sel.volumes_end(); ++iter_sv) + { + Volume* sv = *iter_sv; + // check for each edge whether it has to be copied. + for(size_t i_edge = 0; i_edge < sv->num_edges(); ++i_edge) + { + Edge* e = m_grid.get_edge(sv, i_edge); + + if(m_sel.is_selected(e)) + { + // check the associated vertices through the volumes aaVrtVecVol attachment. + // If at least one has an associated new vertex and if no edge between the + // new vertices already exists, we'll create the new edge. + size_t ind0, ind1; + sv->get_vertex_indices_of_edge(ind0, ind1, i_edge); + Vertex* nv0 = (m_aaVrtVecVol[sv])[ind0]; + Vertex* nv1 = (m_aaVrtVecVol[sv])[ind1]; + + if(nv0 || nv1) + { + // if one vertex has no associated new one, then we use the vertex itself + if(!nv0) + nv0 = sv->vertex(ind0); + if(!nv1) + nv1 = sv->vertex(ind1); + + // create the new edge if it not already exists. + if( ! m_grid.get_edge(nv0, nv1)) + m_grid.create_by_cloning(e, EdgeDescriptor(nv0, nv1), e); + } + } + } + } + + // now we create new faces from selected ones which are connected to + // inner vertices. This allows to preserve old subsets. + // Since we have to make sure that we use the right vertices, + // we have to iterate over the selected volumes and perform all actions on the side-faces + // of those volumes. + + FaceDescriptor fd; + + + for(VolumeIterator iter_sv = m_sel.volumes_begin(); iter_sv != m_sel.volumes_end(); ++iter_sv) + { + Volume* sv = *iter_sv; + // check for each face whether it has to be copied. + for(size_t i_face = 0; i_face < sv->num_faces(); ++i_face) + { + Face* sf = m_grid.get_face(sv, i_face); + + if( m_sel.is_selected(sf)) + { + // check the associated vertices through the volumes aaVrtVecVol attachment. + // If no face between the new vertices already exists, we'll create the new face. + sv->get_vertex_indices_of_face(locVrtInds, i_face); + fd.set_num_vertices(sf->num_vertices()); + + for(size_t i = 0; i < sf->num_vertices(); ++i) + { + Vertex* nVrt = (m_aaVrtVecVol[sv])[locVrtInds[i]]; + if(nVrt) + fd.set_vertex(i, nVrt); + else + fd.set_vertex(i, sv->vertex(locVrtInds[i])); + } + + // if the new face does not already exist, we'll create it + if(!m_grid.get_face(fd)) + m_grid.create_by_cloning(sf, fd, sf); + } + } + } + + // Expand all faces. + // Since volumes are replaced on the fly, we have to take care with the iterator. + // record all new volumes in a vector. This will help to adjust positions later on. + + std::vector newFractureVolumes; + std::vector subsOfNewVolumes; + + VolumeDescriptor vd; + + for(VolumeIterator iter_sv = m_sel.volumes_begin(); iter_sv != m_sel.volumes_end();) + { + Volume* sv = *iter_sv; + ++iter_sv; + + // now expand the fracture faces of sv to volumes. + for(size_t i_side = 0; i_side < sv->num_sides(); ++i_side) + { + // get the local vertex indices of the side of the volume + sv->get_vertex_indices_of_face(locVrtInds, i_side); + + Face* tFace = m_grid.get_side(sv, i_side); + + if(tFace) + { + if(m_aaMarkFaceB[tFace]) + { + Volume* expVol = nullptr; + + if(locVrtInds.size() == 3) + { + size_t iv0 = locVrtInds[0]; + size_t iv1 = locVrtInds[1]; + size_t iv2 = locVrtInds[2]; + + if( ( m_aaVrtVecVol[sv] )[iv0] + && ( m_aaVrtVecVol[sv] )[iv1] + && ( m_aaVrtVecVol[sv] )[iv2] + ) + { + // create a new prism + expVol = *m_grid.create( + PrismDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0), + (m_aaVrtVecVol[sv])[iv2], + (m_aaVrtVecVol[sv])[iv1], + (m_aaVrtVecVol[sv])[iv0])); + } + else if( ( m_aaVrtVecVol[sv] )[iv0] + && ( m_aaVrtVecVol[sv] )[iv1] + ) + { + // create a new Pyramid + expVol = *m_grid.create( + PyramidDescriptor(sv->vertex(iv0), sv->vertex(iv1), + (m_aaVrtVecVol[sv])[iv1], + (m_aaVrtVecVol[sv])[iv0], + sv->vertex(iv2))); + } + else if( ( m_aaVrtVecVol[sv] )[iv1] + && ( m_aaVrtVecVol[sv] )[iv2] + ) + { + // create a new Pyramid + expVol = *m_grid.create( + PyramidDescriptor(sv->vertex(iv1), sv->vertex(iv2), + (m_aaVrtVecVol[sv])[iv2], + (m_aaVrtVecVol[sv])[iv1], + sv->vertex(iv0))); + } + else if( (m_aaVrtVecVol[sv])[iv0] + && (m_aaVrtVecVol[sv])[iv2] + ) + { + // create a new Pyramid + expVol = *m_grid.create( + PyramidDescriptor(sv->vertex(iv2), sv->vertex(iv0), + (m_aaVrtVecVol[sv])[iv0], + (m_aaVrtVecVol[sv])[iv2], + sv->vertex(iv1))); + } + else if( ( m_aaVrtVecVol[sv])[iv0] ) + { + // create a new Tetrahedron + expVol = *m_grid.create( + TetrahedronDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0), + (m_aaVrtVecVol[sv])[iv0])); + } + else if( ( m_aaVrtVecVol[sv])[iv1] ) + { + // create a new Tetrahedron + expVol = *m_grid.create( + TetrahedronDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0), + (m_aaVrtVecVol[sv])[iv1])); + } + else if( ( m_aaVrtVecVol[sv])[iv2] ) + { + // create a new Tetrahedron + expVol = *m_grid.create( + TetrahedronDescriptor(sv->vertex(iv2), sv->vertex(iv1), sv->vertex(iv0), + (m_aaVrtVecVol[sv])[iv2])); + } + else + { + // this code-block should never be entered. If it is entered then + // we either selected the wrong faces (this shouldn't happen), or there + // are selected faces, which have fracture-boundary-vertices only. + // This is the same is if inner fracture edges exists, which are + // connected to two boundary vertices. + // Since we tried to remove those edges above, something went wrong. + // remove the temporary attachments and throw an error +// grid.detach_from_vertices(aVrtVec); +// grid.detach_from_volumes(aVrtVec); +// grid.detach_from_vertices(aAdjMarker); +// grid.detach_from_edges(aAdjMarker); + throw(UGError("Error in ExpandFractures3d. Implementation Error.")); + return false; + } + } + else if ( locVrtInds.size() == 4 ) + { + // newly implemented by Markus to test with Hexahedrons + + size_t iv0 = locVrtInds[0]; + size_t iv1 = locVrtInds[1]; + size_t iv2 = locVrtInds[2]; + size_t iv3 = locVrtInds[3]; + + if( ( m_aaVrtVecVol[sv] )[iv0] + && ( m_aaVrtVecVol[sv] )[iv1] + && ( m_aaVrtVecVol[sv] )[iv2] + && ( m_aaVrtVecVol[sv] )[iv3] + ) + { + // create a new prism + expVol = *m_grid.create( + HexahedronDescriptor( + sv->vertex(iv3), sv->vertex(iv2), + sv->vertex(iv1), sv->vertex(iv0), + (m_aaVrtVecVol[sv])[iv3], + (m_aaVrtVecVol[sv])[iv2], + (m_aaVrtVecVol[sv])[iv1], + (m_aaVrtVecVol[sv])[iv0] + ) + ); + +// m_sh.assign_subset(expVol, m_fracInfosBySubset.at(m_sh.get_subset_index(tFace)).newSubsetIndex); +// +// return true; + } + else if( ( m_aaVrtVecVol[sv] )[iv0] + && ( m_aaVrtVecVol[sv] )[iv1] + + ) + { + // create a new prism + // create a new prism +// expVol = *m_grid.create( +// PrismDescriptor(sv->vertex(iv3),sv->vertex(iv2), sv->vertex(iv1), +// sv->vertex(iv0), +// (m_aaVrtVecVol[sv])[iv1], +// (m_aaVrtVecVol[sv])[iv0]) +// ); + // create a new Prism + /// only used to initialize a prism. for all other tasks you should use VolumeDescripor. + /** + * please be sure to pass the vertices in the correct order: + * v1, v2, v3: bottom-vertices in counterclockwise order (if viewed from the top). + * v4, v5, v6: top-vertices in counterclockwise order (if viewed from the top). + * PrismDescriptor(Vertex* v1, Vertex* v2, Vertex* v3, + Vertex* v4, Vertex* v5, Vertex* v6); + * + */ + expVol = *m_grid.create( + PrismDescriptor( (m_aaVrtVecVol[sv])[iv0], + sv->vertex(iv0), sv->vertex(iv3), + (m_aaVrtVecVol[sv])[iv1], sv->vertex(iv1), sv->vertex(iv2) + ) + ); + + UG_LOG("PRISM 0 1 " << std::endl); + } + else if( ( m_aaVrtVecVol[sv] )[iv0] + && ( m_aaVrtVecVol[sv] )[iv2] + ) + { + UG_LOG("PRISM 0 2 " << std::endl); + expVol = *m_grid.create( + PrismDescriptor( (m_aaVrtVecVol[sv])[iv0], + sv->vertex(iv0), sv->vertex(iv3), + sv->vertex(iv2), (m_aaVrtVecVol[sv])[iv2], sv->vertex(iv1) + ) + ); + +// m_sh.assign_subset(expVol, m_sh.num_subsets()); +// m_sh.assign_subset( sv->vertex(iv0), m_sh.num_subsets()); +// m_sh.assign_subset( sv->vertex(iv1), m_sh.num_subsets()); +// m_sh.assign_subset( sv->vertex(iv2), m_sh.num_subsets()); +// m_sh.assign_subset( sv->vertex(iv3), m_sh.num_subsets()); +// m_sh.assign_subset( (m_aaVrtVecVol[sv])[iv0], m_sh.num_subsets()); +// m_sh.assign_subset( (m_aaVrtVecVol[sv])[iv2], m_sh.num_subsets()); + + + } + else if( ( m_aaVrtVecVol[sv] )[iv0] + && ( m_aaVrtVecVol[sv] )[iv3] + ) + { + UG_LOG("PRISM 0 3 " << std::endl); + + } + else if( ( m_aaVrtVecVol[sv] )[iv1] + && ( m_aaVrtVecVol[sv] )[iv2] + ) + { + + UG_LOG("PRISM 1 2 " << std::endl); + + } + else if( ( m_aaVrtVecVol[sv] )[iv2] + && ( m_aaVrtVecVol[sv] )[iv3] + ) + { + UG_LOG("PRISM 2 3 " << std::endl); + + } + + + + } + else + { + // traditionally only tetrahedrons are supported. This section thus raises an error + // Markus tries to implement also Hexahedra +// grid.detach_from_vertices(aVrtVec); +// grid.detach_from_volumes(aVrtVec); +// grid.detach_from_vertices(aAdjMarker); +// grid.detach_from_edges(aAdjMarker); + throw(UGError("Incomplete implementation error in ExpandFractures3d Arte: Only tetrahedrons are supported in the current implementation, and hexahedra are in development.")); + return false; + } + + if(expVol) + { + + IndexType newSubs = m_fracInfosBySubset.at(m_sh.get_subset_index(tFace)).newSubsetIndex; + + subsOfNewVolumes.push_back( newSubs ); + +// m_sh.assign_subset(expVol, m_fracInfosBySubset.at(m_sh.get_subset_index(tFace)).newSubsetIndex); + m_sh.assign_subset(expVol, newSubs); + + newFractureVolumes.push_back(expVol); + } + } + } + } + + + // now set up a new volume descriptor and replace the volume. + if(vd.num_vertices() != sv->num_vertices()) + vd.set_num_vertices(sv->num_vertices()); + + for(size_t i_vrt = 0; i_vrt < sv->num_vertices(); ++i_vrt) + { + if( (m_aaVrtVecVol[sv])[i_vrt] ) + vd.set_vertex(i_vrt, (m_aaVrtVecVol[sv])[i_vrt]); + else + vd.set_vertex(i_vrt, sv->vertex(i_vrt)); + } + + m_grid.create_by_cloning(sv, vd, sv); + m_grid.erase(sv); + } + + // we have to clean up unused faces and edges. + // note that all selected edges with mark 0 may safley be deleted. - warum? + for(EdgeIterator iter = m_sel.begin(); iter!= m_sel.end();) + { + // take care of the iterator + Edge* e = *iter; + ++iter; + + if( m_aaMarkEdgeVFP[e].getNumberFracEdgesInVertex() == 0 ) + m_grid.erase(e); + } + + // make sure that no unused faces linger around (This should never happen!) + bool foundUnusedFaces = false; + for(FaceIterator iter = m_sel.begin(); iter != m_sel.end();) + { + Face* f = *iter; + ++iter; + + if( ! m_aaMarkFaceB[f] ) + { + foundUnusedFaces = true; + m_grid.erase(f); + } + } + + if(foundUnusedFaces) + { + UG_LOG("WARNING in ExpandFractures3D Arte: Unused faces encountered during cleanup. Removing them...\n"); + } + + if( subsOfNewVolumes.size() != newFractureVolumes.size() ) + { + UG_THROW("andere zahl neue volumes als subdoms " << std::endl); + } + + IndexType nfn = 0; + + for( auto const & nf : newFractureVolumes ) + { + for(size_t iFace = 0; iFace < nf->num_faces(); ++iFace) + { + Face * fac = m_grid.get_face(nf, iFace); + + m_sh.assign_subset( fac, subsOfNewVolumes[nfn] ); + + } + + for(size_t iEdge = 0; iEdge < nf->num_edges(); ++iEdge) + { + Edge* edg = m_grid.get_edge(nf, iEdge); + + m_sh.assign_subset( edg, subsOfNewVolumes[nfn] ); + + } + + for( size_t iVrt = 0; iVrt < nf->num_vertices(); iVrt++ ) + { + Vertex * vrt = nf->vertex(iVrt); + + m_sh.assign_subset( vrt, subsOfNewVolumes[nfn] ); + } + + nfn++; + } + + + + return true; +} + +} /* namespace ug */ diff --git a/ugbase/lib_grid/algorithms/extrusion/ArteExpandFracs3D.h b/ugbase/lib_grid/algorithms/extrusion/ArteExpandFracs3D.h new file mode 100644 index 000000000..385030ad1 --- /dev/null +++ b/ugbase/lib_grid/algorithms/extrusion/ArteExpandFracs3D.h @@ -0,0 +1,302 @@ +/* + * ArteExpandFracs3D.h + * + * Created on: 06.10.2024 + * Author: Markus M. Knodel + * + * * expand fractures using the Arte algorithm, 3D case + * + * Author: Markus Knodel, inspired by Arte from Fuchs and Sebastian Reiters code for fracture expansion without Arte + * + * implementing a class that gives the basic tools for Arte in 3D + * might be templated at a later stage to fulfill 2D and 3D purposes, if suitable + * ( so far 2D case one entire function, not so perfect, but running....) + * + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + +#include +#include +#include +#include "lib_grid/lg_base.h" +#include "expand_layers.h" + //#include "expand_layers_arte.h" + //#include "expand_layers_arte3D.h" +#include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h" +#include "lib_grid/callbacks/callbacks.h" +#include "lib_grid/grid/grid_util.h" +//#include "lib_grid/util/simple_algebra/least_squares_solver.h" + +#include + +#include "support.h" +#include "support3D.h" + + + +#ifndef UGCORE_UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_ARTEEXPANDFRACS3D_H_ +#define UGCORE_UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_ARTEEXPANDFRACS3D_H_ + +namespace ug { + +class ArteExpandFracs3D +{ + +public: + + ArteExpandFracs3D( Grid & grid, SubsetHandler & sh, + std::vector const & fracInfos, + bool useTrianglesInDiamonds, bool establishDiamonds ); + + virtual ~ArteExpandFracs3D(); + + +public: + + bool run(); + +private: + + Grid & m_grid; + SubsetHandler & m_sh; + std::vector m_fracInfos; + bool m_useTrianglesInDiamonds, m_establishDiamonds; + + Grid::VertexAttachmentAccessor m_aaPos; + + // objects for temporary results + FaceDescriptor m_facDescr; + VolumeDescriptor m_volDescr; + +// std::vector m_tmpEdges; // used for temporary results. +// std::vector m_tmpFaces; // used for temporary results. +// std::vector m_tmpVols; // used for temporary results. + + std::vector m_fracInfosBySubset; + + Selector m_sel; + + bool initialize(); + + bool setSelector(); + + bool attachMarkers(); + + bool detachMarkers(); + + using IndexType = unsigned short; + + using AttachedFractFaceEdgeSudo = support::AttachedFractElem; + + using VecAttachedFractFaceEdgeSudo = std::vector; + + using AttachedGenerFaceEdgeSudo = support::AttachedGeneralElem; + + using VecAttachedGenerFaceEdgeSudo = std::vector; + + + using EdgePair = std::pair; + + using VertxFracPropts = support::VertexFracturePropertiesVol; + + using AttVertFracProp = Attachment; + + AttVertFracProp m_aAdjMarkerVFP; + + // TODO FIXME verfehltes Konzept im 3D Fall!!! ERROR + Grid::VertexAttachmentAccessor m_aaMarkVrtVFP; + // TODO FIXME anstatt zu zählen, wieviele fractures angrenzen, was man ja lassen kann, + // vielleicht irgendwo auch gebraucht, muss man vor allem zählen, wieviele subdomains + // von fractures an dem Vertex zusammentreffen!!!!! + // und ob an dem Vertex, wenn er innen ist, eine fracture aufhört, oder weiter geht + // ob also der Vertex "rundum" von fracture faces umgeben ist, oder nur teilweise + // das vertex fracture properties Konzept vom 2D Fall ist also nicht ausreichend + + // AttVertFracProp m_aAdjMarkerVFP; // used to know if an edge is frac edge, suffix vfp misleading.... + + Grid::EdgeAttachmentAccessor m_aaMarkEdgeVFP; + // used to know if an edge is frac edge, suffix vfp misleading.... + + ABool m_aAdjMarkerB; // used to know if an face is frac face + + Grid::FaceAttachmentAccessor m_aaMarkFaceB; + + bool countAndSelectFracBaseNums(); + + std::vector m_originalFractureFaces; + + bool assignOrigFracInfos(); + + +// using AttVrtVec = Attachment >; +// AttVrtVec m_attVrtVec; +// Grid::VolumeAttachmentAccessor m_aaVrtVecVol; + + using AttVecEdge = Attachment>; + using AttVecFace = Attachment>; +// using AttVecVol = Attachment>; + + AttVecEdge m_aAdjInfoEdges; + AttVecFace m_aAdjInfoFaces; +// AttVecVol m_aAdjInfoVols; + + Grid::VertexAttachmentAccessor m_aaVrtInfoAssoEdges; + Grid::VertexAttachmentAccessor m_aaVrtInfoAssoFaces; +// Grid::VertexAttachmentAccessor m_aaVrtInfoAssoVols; + + bool establishNewVrtBase(); + + bool generateVertexInfos(); + + bool loop2EstablishNewVertices(); + +// ug::support::VertexFractureQuadrupel bla(); + +// using VrtxFractrQuadrplArte3D = support::VertexFractureQuadrupel; +// +// using VrtxFractrQuadrplArte3DVec = std::vector; +// +// VrtxFractrQuadrplArte3DVec m_vrtxFractrQuadrplVec; + + using VertFracTrip = support::VertexFractureTripleMF; + + using VecVertFracTrip = std::vector; + + using AttVecVertFracTrip = Attachment; + + AttVecVertFracTrip m_aAdjInfoAVVFT; + + Grid::VertexAttachmentAccessor m_aaVrtInfoFraTri; + + bool checkIfFacesVerticesCoincide( Face * const & facOne, Face * const & facTwo ); + + bool collectFaceVertices( std::vector & facVrt, Face * const & fac ); + + bool createConditionForNewVrtcs(); + + using AttVrtVec = Attachment >; + + AttVrtVec m_attVrtVec; + + Grid::VolumeAttachmentAccessor m_aaVrtVecVol; + + using CrossVertInf = support::CrossingVertexInfoVol; + + std::vector m_vecCrossVrtInf; + + using PairSudoBool = std::pair; + using VecPairSudoBool = std::vector; + + bool isVrtxSurroundedByFracFaces( Vertex * const & vrt, VertxFracPropts & vrtxFracPrps ); +// VecPairSudoBool & sudoSurrounded ); + + // transform to template soon + bool sortElemCircleIsClosed( VecAttachedFractFaceEdgeSudo const & vecAttFac, + VecAttachedFractFaceEdgeSudo & vecSortedFac, + int startFaceIndexUser = -1, +// int endFaceIndexUser = -1, +// IndexType startEdgeIndexUser = -1, +// IndexType endEdgeIndexUser = -1 +// Face * const & startFacUser = nullptr, +// Face * const & endFacUser = nullptr, + Edge * const & startEdgUser = nullptr, + Edge * const & endEdgUser = nullptr + ); + +public: + using VrtxFracProptsStatus = VertxFracPropts::VrtxFracStatus; + +private: +// static_assert< std::is_same< VrtxFracProptsStatus,support::VertexFracturePropertiesVol::VrtxFracStatus>::value ); +// static_assert( std::is_same::value ); +// static_assert( std::is_same::VrtxFracStatus>::value ); + +// template +//// template +//// template +// bool establishNewVertices( Vertex * const & oldVrt ) +// { +// UG_THROW("too general, should not be called presently " << std::endl); +// return false; +// }; + +// template< bool APPLY_GENERAL_SEGMENT_ORDERING, +// ArteExpandFracs3D::VrtxFracProptsStatus vfp, +// std::enable_if::type = true +// > +// bool establishNewVertices( Vertex * const & oldVrt ); + + template< bool APPLY_GENERAL_SEGMENT_ORDERING, + ArteExpandFracs3D::VrtxFracProptsStatus vfp + > + bool establishNewVertices( Vertex * const & oldVrt ); + + +// template< bool APPLY_GENERAL_SEGMENT_ORDERING, +// ArteExpandFracs3D::VrtxFracProptsStatus vfp, +// std::enable_if::type = true +// > +// bool establishNewVertices( Vertex * const & oldVrt ); +// +// template< bool APPLY_GENERAL_SEGMENT_ORDERING, +// ArteExpandFracs3D::VrtxFracProptsStatus vfp, +// typename std::enable_if< std::integral_constant> +// > +// bool establishNewVertices( Vertex * const & oldVrt ); + + bool createNewElements(); + + using AttachedVolumeElemInfo = support::AttachedFullDimElemInfo; + using VecAttachedVolumeElemInfo = std::vector; + using AttVecAttachedVolumeElemInfo = Attachment; + + AttVecAttachedVolumeElemInfo m_aAdjVolElmInfo; + Grid::VertexAttachmentAccessor m_aaVolElmInfo; + + using SegmentVolElmInfo = VecAttachedVolumeElemInfo; + using VecSegmentVolElmInfo = std::vector; + +}; + +// specification has to be declared outside central class context, else compilation error + +template <> +bool ArteExpandFracs3D::establishNewVertices< true, + ArteExpandFracs3D::VrtxFracProptsStatus::oneFracSuDoAtt + >( Vertex * const & oldVrt ); + +template <> +bool ArteExpandFracs3D::establishNewVertices< false, + ArteExpandFracs3D::VrtxFracProptsStatus::oneFracSuDoAtt + >( Vertex * const & oldVrt ); + + +} /* namespace ug */ + +#endif /* UGCORE_UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_ARTEEXPANDFRACS3D_H_ */ diff --git a/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte.cpp b/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte.cpp index c89d0c809..46669c897 100644 --- a/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte.cpp +++ b/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte.cpp @@ -57,7 +57,7 @@ using namespace std; namespace ug{ - +//using namespace ug::support; using VertFracTrip = VertexFractureTriple; @@ -65,7 +65,7 @@ using VertFracTrip = VertexFractureTriple; //using VvftIterator = VecVertFracTrip::iterator; -using AttVrtVec = Attachment >; +using AttVrtVec = Attachment >; //using VertexOfFaceInfo = VertexFractureTriple< std::pair, Face*, std::pair >; // @@ -2533,6 +2533,10 @@ bool ExpandFractures2dArte( Grid& grid, SubsetHandler& sh, vector auto sudoEdg = sh.get_subset_index(*iterEdg); + // sieht quasi doppelt gemoppelt aus + if( sudoEdg != fracInd ) + UG_THROW("Subdomain Index Fehler " << std::endl); + // DEBUG ASSERT TEST static_assert( std::is_same< decltype(sudoEdg), int >::value ); // get vertices of edge, always 2 @@ -2776,7 +2780,7 @@ bool ExpandFractures2dArte( Grid& grid, SubsetHandler& sh, vector for(size_t i_vrt = 0; i_vrt < sf->num_vertices(); ++i_vrt) { - newVrts[i_vrt] = NULL; + newVrts[i_vrt] = nullptr; } // erstmal so tun, als ob keine neuen Vertizes erzeugt werden an den alten Vertizes } @@ -4655,6 +4659,8 @@ bool ExpandFractures2dArte( Grid& grid, SubsetHandler& sh, vector else { +// UG_THROW("in dieser Funktion alt" << std::endl); + // create normal vectors into direction of relevant edges vector3 alongEdgeOne; @@ -4691,7 +4697,7 @@ bool ExpandFractures2dArte( Grid& grid, SubsetHandler& sh, vector } - if( vrtEdgeOneBegin == NULL || vrtEdgeTwoBegin == NULL || vrtEdgeOneEnd == NULL || vrtEdgeTwoEnd == NULL ) + if( vrtEdgeOneBegin == nullptr || vrtEdgeTwoBegin == nullptr || vrtEdgeOneEnd == nullptr || vrtEdgeTwoEnd == nullptr ) { UG_THROW("lauter Nullen vertizes" << std::endl); } @@ -4845,8 +4851,15 @@ bool ExpandFractures2dArte( Grid& grid, SubsetHandler& sh, vector vector3 posNewVrt; - if( cosBetweenNormals > cosinusLim ) +// if( cosBetweenNormals > cosinusLim ) + if( cosBetweenNormals > cosinusLim && numFracsCrossAtVrt < 4 ) +// if( cosBetweenNormals > cosinusLim && subsIndFracOne == subsIndFracTwo ) +// if( false ) { + +// UG_THROW("in dieser Funktion neu cos" << std::endl); + + // dieselben Methoden wie im Fall von einer durchgehenden Kluft an einem Vertex, dort kopieren // bzw Funktion schreiben, die beides macht @@ -4946,6 +4959,9 @@ bool ExpandFractures2dArte( Grid& grid, SubsetHandler& sh, vector } else { + +// UG_THROW("in dieser Funktion neu" << std::endl); + // create normal vectors into direction of relevant edges vector3 alongEdgeOne; @@ -4982,7 +4998,7 @@ bool ExpandFractures2dArte( Grid& grid, SubsetHandler& sh, vector } - if( vrtEdgeOneBegin == NULL || vrtEdgeTwoBegin == NULL || vrtEdgeOneEnd == NULL || vrtEdgeTwoEnd == NULL ) + if( vrtEdgeOneBegin == nullptr || vrtEdgeTwoBegin == nullptr || vrtEdgeOneEnd == nullptr || vrtEdgeTwoEnd == nullptr ) { UG_THROW("lauter Nullen vertizes" << std::endl); } diff --git a/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte3D.cpp b/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte3D.cpp new file mode 100644 index 000000000..78ccf248f --- /dev/null +++ b/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte3D.cpp @@ -0,0 +1,71 @@ +/* + * * expand fractures using the Arte algorithm, 3D case + * + * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt + * Author: Markus Knodel, inspired by Arte from Fuchs and Sebastian Reiters code for fracture expansion without Arte + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + +#include +#include +#include +#include "lib_grid/lg_base.h" +#include "expand_layers.h" +#include "expand_layers_arte.h" +#include "expand_layers_arte3D.h" +#include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h" +#include "lib_grid/callbacks/callbacks.h" +#include "lib_grid/grid/grid_util.h" +//#include "lib_grid/util/simple_algebra/least_squares_solver.h" + +#include + +#include "lib_grid/algorithms/extrusion/ArteExpandFracs3D.h" + +using namespace std; + +namespace ug{ + + +bool ExpandFractures3dArte( Grid& grid, SubsetHandler& sh, + std::vector const & fracInfos, + bool useTrianglesInDiamonds, bool establishDiamonds ) +{ + ArteExpandFracs3D ef3dA ( grid, sh, fracInfos, + useTrianglesInDiamonds, establishDiamonds ); + + return ef3dA.run(); + +} + + + +}// end of namespace + + diff --git a/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte3D.h b/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte3D.h new file mode 100644 index 000000000..a47f93d51 --- /dev/null +++ b/ugbase/lib_grid/algorithms/extrusion/expand_layers_arte3D.h @@ -0,0 +1,53 @@ +/* + * expand fractures using the Arte algorithm, 3D case + * + * Copyright (c) 2011-2015: G-CSC, Goethe University Frankfurt + * Author: Markus Knodel, inspired by Arte from Fuchs and Sebastian Reiters code for fracture expansion without Arte + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + */ + +#ifndef UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_EXPAND_LAYERS_ARTE_3D_H_ +#define UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_EXPAND_LAYERS_ARTE_3D_H_ + +#include +#include "lib_grid/lg_base.h" +#include "expand_layers.h" +#include "expand_layers_arte.h" + +namespace ug +{ + + +bool ExpandFractures3dArte( Grid& grid, SubsetHandler& sh, + std::vector const & fracInfos, + bool useTrianglesInDiamonds, bool establishDiamonds ); + +}// end of namespace + +#endif diff --git a/ugbase/lib_grid/algorithms/extrusion/support.h b/ugbase/lib_grid/algorithms/extrusion/support.h index a641f81ab..a5cfe1ea1 100644 --- a/ugbase/lib_grid/algorithms/extrusion/support.h +++ b/ugbase/lib_grid/algorithms/extrusion/support.h @@ -23,6 +23,12 @@ #include #include +namespace ug +{ + +//namespace support +//{ + ////////////////////////////////////////////////////////////////////////////// @@ -373,7 +379,9 @@ class CrossingVertexInfo ////////////////////////////////////////////////////////////////// +//} +} #endif diff --git a/ugbase/lib_grid/algorithms/extrusion/support3D.h b/ugbase/lib_grid/algorithms/extrusion/support3D.h new file mode 100644 index 000000000..158e038c0 --- /dev/null +++ b/ugbase/lib_grid/algorithms/extrusion/support3D.h @@ -0,0 +1,1343 @@ +/* + * support3D.h + * + * Created on: 31.10.2024 + * Author: Markus Knodel + */ + +#ifndef UGCORE_UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_SUPPORT3D_H_ +#define UGCORE_UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_SUPPORT3D_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "support.h" + +// TODO FIXME +// verschlanken, +// Verdoppelungen weg, +// Namen generalisieren +// nicht vertex spezifisch zB +// sondern lowest dim oder so + +namespace ug +{ + +namespace support +{ + +#if 0 + +template +class ElemInfo +{ +public: + + ElemInfo( ELEMTYP const & elem, INDEX_TYP sudo ) + : m_elem(elem), m_sudo(sudo) + { + } + + + +private: + + ELEMTYP m_elem; + INDEX_TYP m_sudo; + + ElemInfo() {}; + + +}; + +template< +typename FACETYP, +typename NORMALTYP, +typename VOLUMETYP, +typename EDGETYP, +typename INDEX_TYP +> +class VertexFractureQuadrupel +{ +public: + + using ElemInfoFac = ElemInfo; + + //face, normal, volume, edge + +// VertexFractureQuadrupel() +// {}; + + + VertexFractureQuadrupel( ElemInfoFac const & fracFaceInfo, + VOLUMETYP const & attVolume, + NORMALTYP const & normal, + std::pair const & volCutEdges, + std::pair const & volCutEdgeFaces ) + : m_fracFaceInfo(fracFaceInfo), + m_attVolume(attVolume), + m_normal(normal), + m_volCutEdges(volCutEdges), + m_volCutEdgeFaces(volCutEdgeFaces) + { + } + + // todo fixme getter und ggf auch setter, aber vermutlich nur getter implementieren!!! + +private: + +// FACETYP const getFace() const { return m_full; } +// NORMALTYP const getNormal() const { return m_normal; } +// VOLUMETYP const getVolume() const { return m_volume; } +// EDGETYP const getEdge() const { return m_edge; } + + ElemInfo m_fracFaceInfo; + VOLUMETYP m_attVolume; + NORMALTYP m_normal; + std::pair m_volCutEdges; + std::pair m_volCutEdgeFaces; + +//private: +// +// FACETYP m_face; +// NORMALTYP m_normal; +// VOLUMETYP m_volume; +// EDGETYP m_edge; + + VertexFractureQuadrupel() + {}; +}; +#endif + + +template < +typename MANIFELM, +typename LOWDIMELM, +typename INDEX_TXP +> +class AttachedGeneralElem +{ +public: + using PairLowEl = std::pair; + + using AttGenElm = AttachedGeneralElem; + + // for fracture elements + AttachedGeneralElem( MANIFELM const & manifElm, + PairLowEl const & lowElm + ) + : + m_manifElm(manifElm), m_pairLowElm(lowElm) + { + }; + + MANIFELM const getManifElm() const { return m_manifElm;} +// PairLowEl const getLowElm() const { return m_lowElm; } + PairLowEl const getPairLowElm() const { return m_pairLowElm; } + + bool const isNeighboured( AttGenElm const & attElm ) + const + { +// MANIFELM manifElmOther = attElm.getManifElm(); + PairLowEl lowElmOther = attElm.getPairLowElm(); +// INDEX_TXP sudoOther = attElm.getSudo(); + + PairLowEl lowElmThis = this->m_pairLowElm; + + std::vector test; + + test.push_back( lowElmOther.first == lowElmThis.first ); + test.push_back( lowElmOther.second == lowElmThis.first ); + test.push_back( lowElmOther.first == lowElmThis.second ); + test.push_back( lowElmOther.second == lowElmThis.second ); + + INDEX_TXP countCorr = 0; + + for( auto const t : test ) + { + if( t ) + countCorr++; + } + + if( countCorr == 1 ) + return true; + + if( countCorr > 1 ) + UG_THROW("zu viele gleiche Ecken " << std::endl); + + return false; + } + + bool const isNeighbouredAtSpecificSide( AttGenElm const & attElm, + LOWDIMELM const & specificLDE ) + const + { + PairLowEl lowElmOther = attElm.getPairLowElm(); + + PairLowEl lowElmThis = this->m_pairLowElm; + + // test if the specific element is part of at least + // one of the faces + + bool otherFirst = ( lowElmOther.first == specificLDE ); + bool otherSecond = ( lowElmOther.second == specificLDE ); + + bool thisFirst = ( lowElmThis.first == specificLDE ); + bool thisSecond = ( lowElmThis.second == specificLDE ); + + bool isPartOfThisFace = ( thisFirst || thisSecond ); + bool isPartOfOtherFace = ( otherFirst || otherSecond ); + + if( ! isPartOfOtherFace || ! isPartOfThisFace ) + { + UG_LOG("not part of one of the faces " << std::endl); + return false; + } + + if( otherFirst && thisFirst ) + { + if( lowElmOther.first == lowElmThis.first ) + return true; + } + else if( otherFirst && thisSecond ) + { + if( lowElmOther.first == lowElmThis.second ) + return true; + } + else if( otherSecond && thisFirst ) + { + if( lowElmOther.second == lowElmThis.first ) + return true; + } + else if( otherSecond && thisSecond ) + { + if( lowElmOther.second == lowElmThis.second ) + return true; + } + + return false; + } + + bool const testIfEquals( AttGenElm const & attElm ) + const + { + MANIFELM manifElmOther = attElm.getManifElm(); + PairLowEl lowElmOther = attElm.getPairLowElm(); + + if( manifElmOther == this->m_manifElm + && lowElmOther == this->m_pairLowElm + ) + { + return true; + } + + return false; + } + + +protected: + + MANIFELM m_manifElm; + PairLowEl m_pairLowElm; + +}; + + + +/////////////////////////////////////////////////////////////////////////////// + + +// TODO FIXME vertex fracture triplett +// vereinigen mit AttachedGeneralElem !!! davon ableiten!!! +// doppelte Strukturen!!! + + +template < +typename MANIFOLDTYP, // 3D: Face +typename INDEXTYP, // int oder unsinged int oder short oder unsigned short etc +typename FULLDIMTYP, // 3D: Volume +typename SENKRECHTENTYP, // 3D und 2D: ug::vector3 +typename LOWDIMTYP // 3D: Edge (2D nicht benötigt) +> +class VertexFractureTripleMF +: public AttachedGeneralElem +{ + +private: + using AttGenEl = AttachedGeneralElem; + +public: + + using PairLowEl = std::pair; + + VertexFractureTripleMF( MANIFOLDTYP const & manifElm, INDEXTYP sudo, + FULLDIMTYP const & fullElm, + SENKRECHTENTYP const & normal, + PairLowEl const & pairLowElm ) + : //m_manifElm(manifElm), + AttGenEl(manifElm,pairLowElm), + m_sudo(sudo), m_fullElm(fullElm), + m_normal(normal), m_newNormal(normal) +// , +// m_pairLowElm(pairLowElm) + { + }; + +// MANIFOLDTYP const getManifElm() const { return m_manifElm; } + + INDEXTYP const getSudoElm() const { return m_sudo; } + + FULLDIMTYP const getFullElm() const { return m_fullElm; } + + SENKRECHTENTYP const getNormal() const { return m_normal; } + + void setNewNormal( SENKRECHTENTYP const & chNorml ) { m_newNormal = chNorml; } + + SENKRECHTENTYP const getNewNormal() const { return m_newNormal; } + +// PairLowEl const getPairLowElm() const { return m_pairLowElm; } + +private: + +// MANIFOLDTYP m_manifElm; + INDEXTYP m_sudo; + FULLDIMTYP m_fullElm; + SENKRECHTENTYP m_normal; + SENKRECHTENTYP m_newNormal; +// PairLowEl m_pairLowElm; + + VertexFractureTripleMF() + {}; + +}; + + +////////////////////////////////////////////////////////////////// + +// TODO FIXME das muss angepasst werden, ist noch wie für 2D Fall bisher +enum FracTypVol { SingleFrac = 2, TEnd = 3, XCross = 4 }; + +template < typename VRT, typename IndTyp > +class CrossingVertexInfoVol +{ + +public: + + CrossingVertexInfoVol( VRT const & crossVrt, FracTypVol fracTyp ) + : m_crossVrt(crossVrt), m_fracTyp( fracTyp ), + m_vecShiftedVrts(std::vector()) + , m_vecShiftedVrtsWithTypInf(std::vector>()) + , m_numberAtFreeSide(0) + { + } + + VRT getCrossVertex() const { return m_crossVrt; } + + FracTypVol getFracTyp() const { return m_fracTyp; } + + void addShiftVrtx( VRT const & vrt, bool isAtFreeSide = false ) + { + m_vecShiftedVrts.push_back(vrt); + +// if( m_fracTyp == TEnd ) +// { +// std::pair addSVI( vrt, isAtFreeSide ); +// m_vecShiftedVrtsWithTypInf.push_back(addSVI); +// +// if( isAtFreeSide ) +// m_numberAtFreeSide++; +// +// if( m_numberAtFreeSide > 1 ) +// UG_THROW("was ist das fuer ein T Ende" << std::endl); +// } + + } + + void setShiftVrtx( std::vector const & vecVrt ) { m_vecShiftedVrts = vecVrt; } + + std::vector getVecShiftedVrts() const + { + return m_vecShiftedVrts; + } + + std::vector> getVecShiftedVrtsWithTypInfo() const + { +// if( m_fracTyp != TEnd ) +// UG_THROW("fuer Kreuz nicht erlaubt " << std::endl); + + return m_vecShiftedVrtsWithTypInf; + } + + +private: + + VRT m_crossVrt; + std::vector m_vecShiftedVrts; + std::vector> m_vecShiftedVrtsWithTypInf; + FracTypVol m_fracTyp; + IndTyp m_numberAtFreeSide; +}; + +/////////////////////////////////////////////////////////////////////////////// + + + +// info for a vertex: face and attached edges, for fractures only +// derived from the similar class without fracture property! +template < +typename MANIFELM, +typename LOWDIMELM, +typename INDEX_TXP +> +class AttachedFractElem +: public AttachedGeneralElem +// TODO FIXME derive from AttachedGeneralElem +{ +public: + using PairLowEl = std::pair; + + using AttFractElm = AttachedFractElem; + + using AttGenElm = AttachedGeneralElem; + + // for fracture elements + AttachedFractElem( MANIFELM const & manifElm, + PairLowEl & lowElm, + INDEX_TXP sudo ) + : + AttGenElm(manifElm,lowElm), + //m_manifElm(manifElm), m_lowElm(lowElm), + m_sudo(sudo) + { + }; + + +// MANIFELM const getManifElm() const { return m_manifElm;} +// PairLowEl const getLowElm() const { return m_lowElm; } + INDEX_TXP const getSudo() const { return m_sudo; }; + + bool const testIfEquals( AttFractElm const & attElm ) + const + { + bool geomEqu = AttGenElm::testIfEquals(attElm); + +// MANIFELM manifElmOther = attElm.getManifElm(); +// PairLowEl lowElmOther = attElm.getLowElm(); + INDEX_TXP sudoOther = attElm.getSudo(); + +// if( manifElmOther == this->m_manifElm +// && lowElmOther == this->m_lowElm +// && sudoOther == this->m_sudo +// ) +// { +// return true; +// } + + if( geomEqu && sudoOther == this->m_sudo ) + { + return true; + } + + return false; + } + + +// bool const testIfEquals( AttachedFractElem const & attElm ) +// const +// { +// MANIFELM manifElmOther = attElm.getManifElm(); +// PairLowEl lowElmOther = attElm.getLowElm(); +// INDEX_TXP sudoOther = attElm.getSudo(); +// +// if( manifElmOther == this->m_manifElm +// && lowElmOther == this->m_lowElm +// && sudoOther == this->m_sudo +// ) +// { +// return true; +// } +// +// return false; +// } + +// bool const isNeighboured( AttachedFractElem const & attElm ) +// const +// { +//// MANIFELM manifElmOther = attElm.getManifElm(); +// PairLowEl lowElmOther = attElm.getLowElm(); +//// INDEX_TXP sudoOther = attElm.getSudo(); +// +// PairLowEl lowElmThis = this->m_lowElm; +// +// std::vector test; +// +// test.push_back( lowElmOther.first == lowElmThis.first ); +// test.push_back( lowElmOther.second == lowElmThis.first ); +// test.push_back( lowElmOther.first == lowElmThis.second ); +// test.push_back( lowElmOther.second == lowElmThis.second ); +// +// INDEX_TXP countCorr = 0; +// +// for( auto const t : test ) +// { +// if( t ) +// countCorr++; +// } +// +// if( countCorr == 1 ) +// return true; +// +// if( countCorr > 1 ) +// UG_THROW("zu viele gleiche Ecken " << std::endl); +// +// return false; +// } +// +// bool const isNeighbouredAtSpecificSide( AttachedFractElem const & attElm, +// LOWDIMELM const & specificLDE ) +// const +// { +// PairLowEl lowElmOther = attElm.getLowElm(); +// +// PairLowEl lowElmThis = this->m_lowElm; +// +// // test if the specific element is part of at least +// // one of the faces +// +// bool otherFirst = ( lowElmOther.first == specificLDE ); +// bool otherSecond = ( lowElmOther.second == specificLDE ); +// +// bool thisFirst = ( lowElmThis.first == specificLDE ); +// bool thisSecond = ( lowElmThis.second == specificLDE ); +// +// bool isPartOfThisFace = ( thisFirst || thisSecond ); +// bool isPartOfOtherFace = ( otherFirst || otherSecond ); +// +// if( ! isPartOfOtherFace || ! isPartOfThisFace ) +// { +// UG_LOG("not part of one of the faces " << std::endl); +// return false; +// } +// +// if( otherFirst && thisFirst ) +// { +// if( lowElmOther.first == lowElmThis.first ) +// return true; +// } +// else if( otherFirst && thisSecond ) +// { +// if( lowElmOther.first == lowElmThis.second ) +// return true; +// } +// else if( otherSecond && thisFirst ) +// { +// if( lowElmOther.second == lowElmThis.first ) +// return true; +// } +// else if( otherSecond && thisSecond ) +// { +// if( lowElmOther.second == lowElmThis.second ) +// return true; +// } +// +// return false; +// } + +private: +// MANIFELM m_manifElm; +// PairLowEl m_lowElm; + INDEX_TXP m_sudo; + + +}; + +////////////////////////////////////////////////////////////////////////////// + +// class to help count and store a bool and a number of templete type +// comparable to std::pair but more dedicated to the specific aim +// TODO FIXME adapt for 3D case, figure out if inner end, and number of fracs sourrounding +// CAUTION is also used for edges, but still uses +// vertex as indicator - name should be made more flexible + +template< +typename T, +typename ATT_ELEM +> +class VertexFracturePropertiesVol +{ +public: + + using pairTB = std::pair; + using VecPairTB = std::vector; + + +// VertexFracturePropertiesVol( bool isBndFracVertex, T numberCrossingFracsInVertex ) +// : m_isBndFracVertex(isBndFracVertex), m_numberCountedFracsInVertex(numberCrossingFracsInVertex) +// { +// }; + + enum VrtxFracStatus { noFracSuDoAtt = 0, + oneFracSuDoAtt = 1, + twoFracSuDoAtt = 2, + threeFracSuDoAtt = 3 }; + + VertexFracturePropertiesVol() + : m_isBndFracVertex(false), m_numberCountedFracsInVertex(0), + m_status( noFracSuDoAtt ), + m_sudoList( std::vector() ), + m_sudosClosed(VecPairTB()), + m_vecAttElem(std::vector()) +// VertexFracturePropertiesVol( false, 0 ) + { + }; + + void setIsBndFracVertex( bool iBDV = true ) + { + m_isBndFracVertex = iBDV; + } + + void setNumberCrossingFracsInVertex( T const & nCFIV ) + { + m_numberCountedFracsInVertex = nCFIV; + } + + bool getIsBndFracVertex() + { + return m_isBndFracVertex; + } + +// T getCountedNumberFracsInVertex() +// { +// return m_numberCountedFracsInVertex; +// } + + + T getNumberFracEdgesInVertex() + { + return m_numberCountedFracsInVertex; + } + +// T getNumberCrossingFracsInVertex() +// { +// if( m_isBndFracVertex ) +// return m_numberCountedFracsInVertex; +// +// // for inner vertices, each edge passed when +// // fractures are counted along their edges +// // that the vertizes get hit twice for each fracture run +// // only for boundary vertices, this happens only once per fracture +// T multipeInnerHits = 2; +// +// T rest = m_numberCountedFracsInVertex % multipeInnerHits; +// +// if( rest != 0 ) +// { +//// UG_THROW("Expand layers: rest division frac counting not zero " << m_numberCountedFracsInVertex << std::endl); +// +// throw std::runtime_error("error"); +// +// return 0; +// } +// +// return m_numberCountedFracsInVertex / multipeInnerHits; +// } + + VertexFracturePropertiesVol & operator++( int a ) + { + m_numberCountedFracsInVertex++; + return *this; + } + + + VrtxFracStatus getVrtxFracStatus() + { + return m_status; + } + + // caution: returns false if sudo already known, but no problem, feature, not bug + // so never stop the program if false returned here, this is good news + bool addFractSudo( T const & sudo ) + { + bool alreadyInList = false; + + for( auto const & availSudo : m_sudoList ) + { + if( sudo == availSudo ) + alreadyInList = true; + } + + if( ! alreadyInList ) + { + m_sudoList.push_back( sudo ); + } + + return adaptVrtxFracStatus(); + } + + // setter private to avoid abusive use + std::vector const getSudoList() const + { + return m_sudoList; + } + + bool setIsAClosedFracture( T sudoNow, bool isClosedNow ) + { + + T alreadyKnownMult = 0; + + for( auto & suSu : m_sudosClosed ) + { + T & sudoVal = suSu.first; + bool & isClosedVal = suSu.second; + + if( sudoVal == sudoNow ) + { + alreadyKnownMult++; + + UG_LOG("Reassign sudo surround " << std::endl); + + if( isClosedVal != isClosedNow ) + UG_THROW("change property sudo surrounded, why?" << std::endl); + + isClosedVal = isClosedNow; + } + } + + if( alreadyKnownMult == 0 ) + { + pairTB infoSudoSurr( sudoNow, isClosedNow ); + + m_sudosClosed.push_back( infoSudoSurr ); + + } + else if( alreadyKnownMult > 1 ) + { + UG_THROW("zu oft bekannt " << std::endl); + return false; + } + + // check if now correct + + T testKnownFine = 0; + + for( auto const & suSu : m_sudosClosed ) + { + T & sudoVal = suSu.first; + bool & isClosedVal = suSu.second; + + if( sudoVal == sudoNow ) + { + testKnownFine++; + + if( isClosedVal != isClosedNow ) + { + UG_THROW("NOT set property sudo surrounded, why?" << std::endl); + return false; + } + + } + } + + if( testKnownFine == 0 || testKnownFine > 1 ) + { + UG_THROW("immer noch nicht bekannt?" << std::endl); + return false; + + } + + return true; + } + + bool getIsAClosedFracture( T sudoNow ) + { + T foundMultpl = 0; + + bool isClosedReturn = false; + + for( auto const & suSu : m_sudosClosed ) + { + T const & sudoVal = suSu.first; + bool const & isClosedVal = suSu.second; + + if( sudoVal == sudoNow ) + { + foundMultpl++; + isClosedReturn = isClosedVal; + } + } + + if( foundMultpl != 1 ) + { + UG_THROW("not known status closed or not sudo" << std::endl); + return false; + } + + return isClosedReturn; + } + + bool setInfoAllFractureSudosIfClosed( VecPairTB const & sudosClosed ) + { + m_sudosClosed = sudosClosed; + + return true; + } + + VecPairTB const getInfoAllFracSudosIfClosed() const + { + return m_sudosClosed; + } + + // if all open or closed + template + bool const getInfoAllFracturesSameClosedState() const + { + bool allFracsSame = true; + + for( auto const & suSu : m_sudosClosed ) + { + //T const & sudoVal = suSu.first; + bool const & isClosedVal = suSu.second; + + if( isClosedVal != B ) + allFracsSame = false; + } + + return allFracsSame; + } + + bool addAttachedFractElem( ATT_ELEM const & attElem ) + { + bool alreadyKnown = false; + + for( auto const & aE : m_vecAttElem ) + { + if( aE.testIfEquals(attElem) ) + alreadyKnown = true; + } + + if( ! alreadyKnown ) + m_vecAttElem.push_back(attElem); + + // returns true if ads it, false if no need as known + return ! alreadyKnown; + } + + std::vector const & getAllAttachedFractElems() + const + { + return m_vecAttElem; + } + +private: + bool m_isBndFracVertex; + T m_numberCountedFracsInVertex; + + VrtxFracStatus m_status; + + static VrtxFracStatus constexpr m_maxStatus = VrtxFracStatus::threeFracSuDoAtt; + + std::vector m_sudoList; + + // better private, to avoid confusion + bool setVrtxFracStatus( VrtxFracStatus status ) + { + m_status = status; + + if( status < noFracSuDoAtt || status > threeFracSuDoAtt ) + return false; + + return true; + } + + bool adaptVrtxFracStatus() + { + auto sudosNum = m_sudoList.size(); + if( sudosNum > static_cast( m_maxStatus ) ) + { + UG_THROW("zu viele subdomains crossing in one Punkt" << std::endl); + return false; + } + + m_status = static_cast(sudosNum); + + return true; + } + + VecPairTB m_sudosClosed; + + bool setSudoList( std::vector const & sudoList ) + { + m_sudoList = sudoList; + + return true; + } + + std::vector m_vecAttElem; + +}; + + +// intention, explained for volume: +template +class AttachedFullDimElemInfo +{ + +public: + + using AttachedFractManifElemInfo = AttachedFractElem; + using AttachedGenerManifElemInfo = AttachedGeneralElem; + + using VecAttachedFractManifElemInfo = std::vector; + using VecAttachedGenerManifElemInfo = std::vector; + + using AttFullDimElmInfo = AttachedFullDimElemInfo; + + AttachedFullDimElemInfo( FULLDIM_ELEM const & fullDimElm ) + : m_fullDimElm(fullDimElm), + m_elementMarked(false), + m_vecFractManifElm(VecAttachedFractManifElemInfo()), +// m_vecFractManifElmTouchInfo(VecAttFractManifElmTouchInf()), +// m_allSidesTouched(false), + m_vecGenerManifElm(VecAttachedGenerManifElemInfo()) +// m_vecGenerManifElmTouchInfo(VecAttFractManifElmTouchInf()) +// m_vecInnerSegmentManifElm(VecAttachedGenerManifElemInfo()) + { + } + + FULLDIM_ELEM const getFulldimElem() const + { + return m_fullDimElm; + } + + bool const hasSameFulldimElem( AttFullDimElmInfo const & otherFullDimElmInf ) + const + { + FULLDIM_ELEM otherFullDimElm = otherFullDimElmInf.getFulldimElem(); + + return ( otherFullDimElm == m_fullDimElm ); + } + + bool const isMarked() const { return m_elementMarked; } + + void markIt() { m_elementMarked = true; } // to allow in loop to be start element + + bool const hasFracture() const { return ( m_vecFractManifElm.size() > 0 ); } + + bool addFractManifElem( AttachedFractManifElemInfo const & manifFractElm, Grid & grid ) + { + return addManifElem( manifFractElm, m_vecFractManifElm, grid ); + } + + bool addGenerManifElem( AttachedGenerManifElemInfo const & manifGenerElm, Grid & grid ) + { +// static_assert(std::is_same::value); + + return addManifElem( manifGenerElm, m_vecGenerManifElm, grid ); + } + + // necessary to avoid stupid casting from derived class AttachedFractManifElemInfo + // else, addGenerManifElem would also eat objects of derived class + // however not it only accepts explicit base class objects + template + bool addGenerManifElem( NOGEN const & noGener, Grid & grid ) + = delete; + + +// // will form new "surface" of next inner round in a segment +// bool addInnerSegmentManifElem( AttachedGenerManifElemInfo const & manifInnerSegmentElm, Grid & grid ) +// { +// return addManifElem( manifInnerSegmentElm, m_vecInnerSegmentManifElm, grid ); +// } +// +// template +// bool addInnerSegmentElem( NOGEN const & noGener, Grid & grid ) +// = delete; + + + +// bool addFractureManifElem( AttachedFractManifElemInfo const & manifElm, Grid & grid ) +// { +// // Caution: first test if manifold elem is at all part of the fulldim elem manifols +// // if not, return false directly +// +// if( ! fullDimElmContainsManif( manifElm.getManifElm(), grid ) ) +// return false; +// +// // if manif elem is in principle part of the fulldim elem manifolds, +// // then we need to check if it is already integrated +// +// bool hasElemAlready = false; +// +// +// for( auto const & me : m_vecFractManifElm ) +// { +// if( manifElm.testIfEquals(me) ) +// { +// hasElemAlready = true; +// break; +// } +// } +// +// if( ! hasElemAlready ) +// { +// m_vecFractManifElm.push_back( manifElm ); +// +//// AttFractManifElmTouchInf pamei( manifElm, false ); +//// +//// m_vecFractManifElmTouchInfo.push_back(pamei); +// } +// +// return ! hasElemAlready; +// } + + VecAttachedFractManifElemInfo const getVecFractManifElem() const + { + return m_vecFractManifElm; + } + + VecAttachedGenerManifElemInfo const getVecGenerManifElem() const + { + return m_vecGenerManifElm; + } + + bool const searchGenerManifElem( AttachedGenerManifElemInfo const & manifGenerElemOther, bool eraseFound = true ) + { + bool found = searchManifElem( manifGenerElemOther, m_vecGenerManifElm, eraseFound ); + + if( found && eraseFound ) + { + m_elementMarked = true; + } + + return found; + } + + bool const testFullDimElmNeighbour( AttFullDimElmInfo const & attFullDimElmInfOther, bool eraseFoundManif = true ) + { + VecAttachedGenerManifElemInfo const & vecGenerManifElmOther = attFullDimElmInfOther.getVecGenerManifElem(); + + bool manifNeighbored = false; + + for( AttachedGenerManifElemInfo const & generManifElemOther : vecGenerManifElmOther ) + { + if( searchManifElem( generManifElemOther, m_vecGenerManifElm, eraseFoundManif ) ) + manifNeighbored = true; + + } + + if( manifNeighbored && eraseFoundManif ) + { + m_elementMarked = true; + } + + return manifNeighbored; + } + + + template + bool searchGenerManifElem( NOGEN const & manifGenerElemOther, bool eraseFound ) = delete; + +// bool const searchFractManifElem( AttachedFractManifElemInfo const & manifFractElemOther, bool eraseFound = true ) + bool const searchFractManifElem( AttachedFractManifElemInfo const & manifFractElemOther, bool shiftToGeneral = true ) + { + bool found = searchManifElem( manifFractElemOther, m_vecFractManifElm, shiftToGeneral ); + + if( found && shiftToGeneral ) + { + // for the case that a fracture is not closed at a vertex, shift the in principle + // fracture vertex to the gerneral vertices, as useless for segmente construction + // and useless for expansion + + MANIFELM const & manifel = manifFractElemOther.getManifElm(); + typename AttachedGenerManifElemInfo::PairLowEl const & pairLowEl = manifFractElemOther.getPairLowElm(); + + AttachedGenerManifElemInfo agmei( manifel, pairLowEl ); + + m_vecGenerManifElm.push_back(agmei); + } + + return found; + + } + +// bool const searchInnerSegmentManifElem( AttachedGenerManifElemInfo const & manifInnerSegmElemOther, bool eraseFound = true ) +// { +// return searchManifElem( manifInnerSegmElemOther, m_vecInnerSegmentManifElm, eraseFound ); +// } + +// bool const tryToTouchManifElem( AttachedFractManifElemInfo const & manifElemOther ) const +// { +//// bool managed2Touch = false; +// +// for( typename VecAttachedFractManifElemInfo::iterator afeiIter = m_vecFractManifElm.begin(); +// afeiIter != m_vecFractManifElm.end(); +// afeiIter++ +// ) +// { +// AttachedFractManifElemInfo & manifElmTest = *afeiIter; +// +// if( manifElemOther.testIfEquals(manifElmTest) ) +// { +// managed2Touch = true; +// +// m_vecFractManifElm.erase(afeiIter); +// +// return true; +// } +// } +// +// return false; +// +//// return managed2Touch; +// +//// for( auto & ameti : m_vecFractManifElmTouchInfo ) +//// { +//// AttachedFractManifElemInfo & manifElmTest = ameti.first; +//// bool & alreadyTouched = ameti.second; +//// +//// if( ! alreadyTouched ) +//// { +//// if( manifElemOther.testIfEquals( manifElmTest ) ) +//// { +//// alreadyTouched = true; +//// } +//// } +//// +//// if( alreadyTouched ) +//// { +//// managed2Touch = true; +//// break; +//// } +//// } +// +// return managed2Touch; +// } + +// VecAttachedFractManifElemInfo const getAlreadyTouchedManifElems() const +// { +// VecAttachedFractManifElemInfo alreadyTouchedManifElms; +// +// for( const auto & ameti : m_vecFractManifElmTouchInfo ) +// { +// if( ameti.second ) +// alreadyTouchedManifElms.push_back( ameti ); +// } +// +// return alreadyTouchedManifElms; +// } + +// VecAttachedFractManifElemInfo const getSoFarUnTouchedManifElems() const +// { +// VecAttachedFractManifElemInfo unTouchedManifElms; +// +// for( const auto & ameti : m_vecFractManifElmTouchInfo ) +// { +// if( ! ameti.second ) +// unTouchedManifElms.push_back( ameti ); +// } +// +// return unTouchedManifElms; +// } + +// bool const testIfAllSidesTouched() const +// { +// if( m_allSidesTouched ) +// return true; +// +// bool allSidesTouched = true; +// +// for( const auto & ameti : m_vecFractManifElmTouchInfo ) +// { +// if( ! ameti.second ) +// { +// allSidesTouched = false; +// } +// } +// +// m_allSidesTouched = allSidesTouched; +// +// return m_allSidesTouched; +// } + +private: + +// bool m_allSidesTouched; + + FULLDIM_ELEM m_fullDimElm; + + bool m_elementMarked; + + VecAttachedFractManifElemInfo m_vecFractManifElm; + +// using AttFractManifElmTouchInf = std::pair; +// using VecAttFractManifElmTouchInf = std::vector; +// +// VecAttFractManifElmTouchInf m_vecFractManifElmTouchInfo; + + VecAttachedGenerManifElemInfo m_vecGenerManifElm; + +// VecAttachedGenerManifElemInfo m_vecInnerSegmentManifElm; + +// using AttGenerManifElmTouchInf = std::pair; +// using VecAttGenerManifElmTouchInf = std::vector; +// +// VecAttFractManifElmTouchInf m_vecGenerManifElmTouchInfo; + + +// template< +//// typename std::enable_if::value, true>, +//// typename std::enable_if::value, true> +//// typename std::enable_if::value, true>, +// typename std::enable_if::value,bool> +// > + //std::enable_if::value> + template + < + typename = std::enable_if::value>, + typename = std::enable_if::value> +// typename = std::enable_if::value>, +// typename = std::enable_if::value> + > + bool fullDimElmContainsManif( MANIFELM const & manifEl, Grid & grid ) + { +// bool contained = false; + + for( INDEX_TXP iFac = 0; iFac < m_fullDimElm->num_faces(); iFac++ ) + { + +// static_assert(std::is_same::value); + + Face * fac = grid.get_face(m_fullDimElm,iFac); + + if( fac == manifEl ) + { + return true; +// contained = true; +// break; + } + } + +// return contained; + return false; + } + + + template + bool const searchManifElem( ATT_MANIF_ELM_INF const & manifElemOther, + std::vector & memVecManifElem, + bool eraseFound = true ) + const + { + + for( typename std::vector::iterator afeiIter = memVecManifElem.begin(); + afeiIter != memVecManifElem.end(); + afeiIter++ + ) + { + ATT_MANIF_ELM_INF & manifElmTest = *afeiIter; + + if( manifElemOther.testIfEquals(manifElmTest) ) + { + + if( eraseFound ) + memVecManifElem.erase(afeiIter); + + return true; + } + } + + return false; + } + + + + template + bool addManifElem( ATT_MANIF_ELM_INFO const & manifElm, + std::vector & memVecManifElm, + Grid & grid ) + { + // Caution: first test if manifold elem is at all part of the fulldim elem manifols + // if not, return false directly + + if( ! fullDimElmContainsManif( manifElm.getManifElm(), grid ) ) + return false; + + // if manif elem is in principle part of the fulldim elem manifolds, + // then we need to check if it is already integrated + + for( auto const & me : memVecManifElm ) + { + if( manifElm.testIfEquals(me) ) + { + return false; + } + } + + // not contained so far, but part of the manifolds of the fulldim elem + memVecManifElm.push_back(manifElm); + + return true; + + } + + +}; + +#if 0 +template +bool switchFulldimInfo( VEC_AVEI & vecAttVolElemInfoCop, + VEC_AVEI const & vecAttVolElemInfo, + VEC_AVEI & segmentAVEI, + OPERATION opera, + INDX_TYP switchInd = 0 + ) +{ + auto & startVolInfoThisSegment = vecAttVolElemInfoCop[switchInd]; + + auto const & startVol = startVolInfoThisSegment.opera(); + + for( auto & possibleOrigVolInfo : vecAttVolElemInfo ) + { + auto const & possVol = possibleOrigVolInfo.opera(); + + if( possVol == startVol ) + { + segmentAVEI().push_back(possibleOrigVolInfo); + break; + } + } + + if( segmentAVEI().size() != 1 ) + { + UG_LOG("No start volume reconstructible " << std::endl); + UG_THROW("No start volume reconstructible " << std::endl); + return false; + } + + if( ! vecAttVolElemInfoCop.erase( vecAttVolElemInfoCop.begin() + switchInd ) ) + return false; + + return true; + +} +#endif + +} + +} + +#endif /* UGCORE_UGBASE_LIB_GRID_ALGORITHMS_EXTRUSION_SUPPORT3D_H_ */