From aaf09124917f5e54a80b9a19c8a9345b5233c9d0 Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 17 Dec 2008 15:56:37 +0000 Subject: [PATCH] MEDMEM Industrialization 2008 calculate barycentre of intersection --- src/INTERP_KERNEL/IntersectorTetra.hxx | 26 +++- src/INTERP_KERNEL/IntersectorTetra.txx | 41 +++++- .../TriangulationIntersector.hxx | 21 +++- .../TriangulationIntersector.txx | 118 +++++++++++------- 4 files changed, 149 insertions(+), 57 deletions(-) diff --git a/src/INTERP_KERNEL/IntersectorTetra.hxx b/src/INTERP_KERNEL/IntersectorTetra.hxx index a419c43e6..5b4dc72c7 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.hxx +++ b/src/INTERP_KERNEL/IntersectorTetra.hxx @@ -1,3 +1,21 @@ +// Copyright (C) 2007-2008 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// #ifndef __INTERSECTORTETRA_HXX__ #define __INTERSECTORTETRA_HXX__ @@ -175,8 +193,8 @@ namespace INTERP_KERNEL /** * \brief Class calculating the volume of intersection between a tetrahedral target element and - * source elements with triangular or quadratilateral faces. - * + * source elements with triangular or quadratilateral faces. Optionally it computes barycentre + * of intersection */ template class INTERPKERNEL_EXPORT IntersectorTetra : public TargetIntersector @@ -191,7 +209,7 @@ namespace INTERP_KERNEL ~IntersectorTetra(); - double intersectSourceCell(typename MyMeshType::MyConnType srcCell); + double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0); private: @@ -224,7 +242,7 @@ namespace INTERP_KERNEL /// reference to the source mesh const MyMeshType& _srcMesh; - + }; /** diff --git a/src/INTERP_KERNEL/IntersectorTetra.txx b/src/INTERP_KERNEL/IntersectorTetra.txx index b501cec45..4ea8daf51 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.txx +++ b/src/INTERP_KERNEL/IntersectorTetra.txx @@ -1,3 +1,21 @@ +// Copyright (C) 2007-2008 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// #ifndef __INTERSECTORTETRA_TXX__ #define __INTERSECTORTETRA_TXX__ @@ -7,8 +25,9 @@ #include "TransformedTriangle.hxx" #include "MeshUtils.hxx" #include "VectorUtils.hxx" -#include "CellModel.hxx" +#include "UniversalCellModel.hxx" #include "Log.hxx" +#include "UnitTetraIntersectionBary.hxx" #include #include @@ -83,11 +102,12 @@ namespace INTERP_KERNEL * with triangulated faces to avoid having to recalculate these. * * @param element global number of the source element (1 <= srcCell < # source cells) + * @param element bary centre of cells intersection */ template - double IntersectorTetra::intersectSourceCell(typename MyMeshType::MyConnType element) + double IntersectorTetra::intersectSourceCell(typename MyMeshType::MyConnType element, + double* baryCentre) { - //{ could be done on outside? // check if we have planar tetra element if(_t->determinant() == 0.0) @@ -98,8 +118,9 @@ namespace INTERP_KERNEL } // get type of cell - NormalizedCellType normCellType=_srcMesh.getTypeOfElement(element); - const CellModel& cellModelCell=CellModel::getCellModel(normCellType); +// NormalizedCellType normCellType=_srcMesh.getTypeOfElement(element); +// const CellModel& cellModelCell=CellModel::getCellModel(normCellType); + UniversalCellModel cellModelCell( _srcMesh, element ); unsigned nbOfNodes4Type=cellModelCell.getNumberOfNodes(); // halfspace filtering bool isOutside[8] = {true, true, true, true, true, true, true, true}; @@ -134,6 +155,9 @@ namespace INTERP_KERNEL if(!isTargetOutside) { + /// calculator of intersection barycentre + UnitTetraIntersectionBary baryCalculator; + for(unsigned ii = 0 ; ii < cellModelCell.getNumberOfSons() ; ++ii) { NormalizedCellType faceType = cellModelCell.getSonType(ii); @@ -154,6 +178,8 @@ namespace INTERP_KERNEL TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); calculateVolume(tri, key); totalVolume += _volumes[key]; + if ( baryCentre ) + baryCalculator.addSide( tri ); } else { // count negative as face has reversed orientation totalVolume -= _volumes[key]; @@ -212,6 +238,11 @@ namespace INTERP_KERNEL } delete [] faceNodes; } + + if ( baryCentre ) { + baryCalculator.getBary( baryCentre ); + _t->reverseApply( baryCentre, baryCentre ); + } } delete [] cellNodes; // reset if it is very small to keep the matrix sparse diff --git a/src/INTERP_KERNEL/TriangulationIntersector.hxx b/src/INTERP_KERNEL/TriangulationIntersector.hxx index 04a5ed65c..56704a045 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.hxx +++ b/src/INTERP_KERNEL/TriangulationIntersector.hxx @@ -1,3 +1,21 @@ +// Copyright (C) 2007-2008 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// #ifndef _TRIANGULATION_INTERSECTOR_HXX_ #define _TRIANGULATION_INTERSECTOR_HXX_ @@ -16,7 +34,8 @@ namespace INTERP_KERNEL public: TriangulationIntersector(const MyMeshType& mesh_A, const MyMeshType& mesh_B, double dimCaracteristic, double precision, double medianPlane, int printLevel); - double intersectCells(ConnType icell_A, ConnType icell_B, int nb_NodesA, int nb_NodesB); + double intersectCells(ConnType icell_A, ConnType icell_B, int nb_NodesA, int nb_NodesB, + double* baryCentre=0); private : const ConnType *_connectA; const ConnType *_connectB; diff --git a/src/INTERP_KERNEL/TriangulationIntersector.txx b/src/INTERP_KERNEL/TriangulationIntersector.txx index ac4114189..9a2ae9653 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.txx +++ b/src/INTERP_KERNEL/TriangulationIntersector.txx @@ -1,3 +1,21 @@ +// Copyright (C) 2007-2008 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// #ifndef __TRIANGULATIONINTERSECTOR_TXX__ #define __TRIANGULATIONINTERSECTOR_TXX__ @@ -31,72 +49,78 @@ namespace INTERP_KERNEL template double TriangulationIntersector::intersectCells(ConnType icell_A, ConnType icell_B, - int nb_NodesA, int nb_NodesB) + int nb_NodesA, int nb_NodesB, + double* baryCentre) { double result = 0.; - int orientation = 1; - + int orientation = 1; + //Obtain the coordinates of A and B std::vector Coords_A (SPACEDIM*nb_NodesA); std::vector Coords_B (SPACEDIM*nb_NodesB); for (int idim=0; idim::coo2C(_connectA[OTT::conn2C(_connIndexA[OTT::ind2C(icell_A)]+i_A)])+idim]; - for (ConnType i_B=0; i_B::coo2C(_connectB[OTT::conn2C(_connIndexB[OTT::ind2C(icell_B)]+i_B)])+idim]; - } - + { + for (ConnType i_A=0; i_A::coo2C(_connectA[OTT::conn2C(_connIndexA[OTT::ind2C(icell_A)]+i_A)])+idim]; + for (ConnType i_B=0; i_B::coo2C(_connectB[OTT::conn2C(_connIndexB[OTT::ind2C(icell_B)]+i_B)])+idim]; + } + //project cells A and B on the median plane and rotate the median plane if(SPACEDIM==3) orientation = projection(Coords_A, Coords_B, nb_NodesA, nb_NodesB, - PlanarIntersector::_dimCaracteristic * PlanarIntersector::_precision, - PlanarIntersector::_medianPlane, PlanarIntersector::_doRotate); - + PlanarIntersector::_dimCaracteristic * PlanarIntersector::_precision, + PlanarIntersector::_medianPlane, PlanarIntersector::_doRotate); + //DEBUG PRINTS if(PlanarIntersector::_printLevel >= 3) - { - std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl; - std::cout << std::endl << "icell_A= " << icell_A << ", nb nodes A= " << nb_NodesA << std::endl; - for(int i_A =0; i_A< nb_NodesA; i_A++) - {for (int idim=0; idim inter; - INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[SPACEDIM*i_A],&Coords_A[SPACEDIM*(i_A+1)], - &Coords_B[0],&Coords_B[SPACEDIM*i_B],&Coords_B[SPACEDIM*(i_B+1)], - inter, PlanarIntersector::_dimCaracteristic, - PlanarIntersector::_precision); - ConnType nb_inter=((ConnType)inter.size())/2; - if(nb_inter >3) inter=reconstruct_polygon(inter); - for(ConnType i = 1; i(&inter[0],&inter[2*i],&inter[2*(i+1)],area); - result +=0.5*fabs(area[0]); - } - //DEBUG prints - if(PlanarIntersector::_printLevel >= 3) - { - std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl; - for(ConnType i=0; i< nb_inter; i++) - {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; } - } - } + std::vector inter; + INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[SPACEDIM*i_A],&Coords_A[SPACEDIM*(i_A+1)], + &Coords_B[0],&Coords_B[SPACEDIM*i_B],&Coords_B[SPACEDIM*(i_B+1)], + inter, PlanarIntersector::_dimCaracteristic, + PlanarIntersector::_precision); + ConnType nb_inter=((ConnType)inter.size())/2; + if(nb_inter >3) inter=reconstruct_polygon(inter); + for(ConnType i = 1; i(&inter[0],&inter[2*i],&inter[2*(i+1)],area); + result +=0.5*fabs(area[0]); + } + if ( baryCentre ) { + std::vector Bary=INTERP_KERNEL::bary_poly(inter); + baryCentre[0] = Bary[0]; + baryCentre[1] = Bary[1]; + } + //DEBUG prints + if(PlanarIntersector::_printLevel >= 3) + { + std::cout << std::endl << "Number of nodes of the intersection = "<< nb_inter << std::endl; + for(ConnType i=0; i< nb_inter; i++) + {for (int idim=0; idim<2; idim++) std::cout << inter[2*i+idim] << " "; std::cout << std::endl; } + } } - + } + //DEBUG PRINTS if(PlanarIntersector::_printLevel >= 3) std::cout << std::endl <<"Intersection area = " << result << std::endl; - + return orientation*result; } } -- 2.39.2