+// 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__
/**
* \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 MyMeshType>
class INTERPKERNEL_EXPORT IntersectorTetra : public TargetIntersector<typename MyMeshType::MyConnType>
~IntersectorTetra();
- double intersectSourceCell(typename MyMeshType::MyConnType srcCell);
+ double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
private:
/// reference to the source mesh
const MyMeshType& _srcMesh;
-
+
};
/**
+// 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__
#include "TransformedTriangle.hxx"
#include "MeshUtils.hxx"
#include "VectorUtils.hxx"
-#include "CellModel.hxx"
+#include "UniversalCellModel.hxx"
#include "Log.hxx"
+#include "UnitTetraIntersectionBary.hxx"
#include <cmath>
#include <cassert>
* 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<class MyMeshType>
- double IntersectorTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element)
+ double IntersectorTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element,
+ double* baryCentre)
{
-
//{ could be done on outside?
// check if we have planar tetra element
if(_t->determinant() == 0.0)
}
// 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<MyMeshType> cellModelCell( _srcMesh, element );
unsigned nbOfNodes4Type=cellModelCell.getNumberOfNodes();
// halfspace filtering
bool isOutside[8] = {true, true, true, true, true, true, true, true};
if(!isTargetOutside)
{
+ /// calculator of intersection barycentre
+ UnitTetraIntersectionBary baryCalculator;
+
for(unsigned ii = 0 ; ii < cellModelCell.getNumberOfSons() ; ++ii)
{
NormalizedCellType faceType = cellModelCell.getSonType(ii);
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];
}
delete [] faceNodes;
}
+
+ if ( baryCentre ) {
+ baryCalculator.getBary( baryCentre );
+ _t->reverseApply( baryCentre, baryCentre );
+ }
}
delete [] cellNodes;
// reset if it is very small to keep the matrix sparse
+// 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_
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;
+// 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__
template<class MyMeshType>
double TriangulationIntersector<MyMeshType>::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<double> Coords_A (SPACEDIM*nb_NodesA);
std::vector<double> Coords_B (SPACEDIM*nb_NodesB);
for (int idim=0; idim<SPACEDIM; idim++)
- {
- for (ConnType i_A=0; i_A<nb_NodesA; i_A++)
- Coords_A[SPACEDIM*i_A+idim] = _coordsA[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectA[OTT<ConnType,numPol>::conn2C(_connIndexA[OTT<ConnType,numPol>::ind2C(icell_A)]+i_A)])+idim];
- for (ConnType i_B=0; i_B<nb_NodesB; i_B++)
- Coords_B[SPACEDIM*i_B+idim] = _coordsB[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectB[OTT<ConnType,numPol>::conn2C(_connIndexB[OTT<ConnType,numPol>::ind2C(icell_B)]+i_B)])+idim];
- }
-
+ {
+ for (ConnType i_A=0; i_A<nb_NodesA; i_A++)
+ Coords_A[SPACEDIM*i_A+idim] = _coordsA[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectA[OTT<ConnType,numPol>::conn2C(_connIndexA[OTT<ConnType,numPol>::ind2C(icell_A)]+i_A)])+idim];
+ for (ConnType i_B=0; i_B<nb_NodesB; i_B++)
+ Coords_B[SPACEDIM*i_B+idim] = _coordsB[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectB[OTT<ConnType,numPol>::conn2C(_connIndexB[OTT<ConnType,numPol>::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<MyMeshType>::_dimCaracteristic * PlanarIntersector<MyMeshType>::_precision,
- PlanarIntersector<MyMeshType>::_medianPlane, PlanarIntersector<MyMeshType>::_doRotate);
-
+ PlanarIntersector<MyMeshType>::_dimCaracteristic * PlanarIntersector<MyMeshType>::_precision,
+ PlanarIntersector<MyMeshType>::_medianPlane, PlanarIntersector<MyMeshType>::_doRotate);
+
//DEBUG PRINTS
if(PlanarIntersector<MyMeshType>::_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<SPACEDIM; idim++) std::cout << Coords_A[SPACEDIM*i_A+idim] << " "; std::cout << std::endl;}
- std::cout << std::endl << "icell_B= " << icell_B << ", nb nodes B= " << nb_NodesB << std::endl;
- for(int i_B =0; i_B< nb_NodesB; i_B++)
- {for (int idim=0; idim<SPACEDIM; idim++) std::cout << Coords_B[SPACEDIM*i_B+idim]<< " "; std::cout << std::endl; }
- }
-
+ {
+ 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<SPACEDIM; idim++) std::cout << Coords_A[SPACEDIM*i_A+idim] << " "; std::cout << std::endl;}
+ std::cout << std::endl << "icell_B= " << icell_B << ", nb nodes B= " << nb_NodesB << std::endl;
+ for(int i_B =0; i_B< nb_NodesB; i_B++)
+ {for (int idim=0; idim<SPACEDIM; idim++) std::cout << Coords_B[SPACEDIM*i_B+idim]<< " "; std::cout << std::endl; }
+ }
+
//Compute the intersection area
double area[SPACEDIM];
for(ConnType i_A = 1; i_A<nb_NodesA-1; i_A++)
+ {
+ for(ConnType i_B = 1; i_B<nb_NodesB-1; i_B++)
{
- for(ConnType i_B = 1; i_B<nb_NodesB-1; i_B++)
- {
- std::vector<double> 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<MyMeshType>::_dimCaracteristic,
- PlanarIntersector<MyMeshType>::_precision);
- ConnType nb_inter=((ConnType)inter.size())/2;
- if(nb_inter >3) inter=reconstruct_polygon(inter);
- for(ConnType i = 1; i<nb_inter-1; i++)
- {
- INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
- result +=0.5*fabs(area[0]);
- }
- //DEBUG prints
- if(PlanarIntersector<MyMeshType>::_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<double> 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<MyMeshType>::_dimCaracteristic,
+ PlanarIntersector<MyMeshType>::_precision);
+ ConnType nb_inter=((ConnType)inter.size())/2;
+ if(nb_inter >3) inter=reconstruct_polygon(inter);
+ for(ConnType i = 1; i<nb_inter-1; i++)
+ {
+ INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
+ result +=0.5*fabs(area[0]);
+ }
+ if ( baryCentre ) {
+ std::vector<double> Bary=INTERP_KERNEL::bary_poly(inter);
+ baryCentre[0] = Bary[0];
+ baryCentre[1] = Bary[1];
+ }
+ //DEBUG prints
+ if(PlanarIntersector<MyMeshType>::_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<MyMeshType>::_printLevel >= 3)
std::cout << std::endl <<"Intersection area = " << result << std::endl;
-
+
return orientation*result;
}
}
