[Dune] [Dune-Commit] dune-grid r6536 - trunk/dune/grid/uggrid
Martin Nolte
nolte at mathematik.uni-freiburg.de
Sat Mar 27 08:03:04 CET 2010
Hi Oli,
first of all thanks for the additional documentation. I just wanted to
point out that all coordinates have to be set for the inside element.
The outernormal is computed as the (inverse) Piola transformation of the
normal to the reference element, i.e.,
JInvT = jacobianInverseTransposed( x )
det = integrationElement( x )
normal = (1/det) * JInvT * refNormal( face )
In order to compute jacobianInverseTransposed and integrationElement,
you need all corners.
In the case of conforming intersections, you might also be interested in
another feature of the GenericGeometry implementation: It is possible to
obtain traces, i.e., to restrict a geometry to a subentity. Though the
current implementation just creates a new geometry, this could be sped
up in the future.
Yours,
Martin
On 03/26/2010 10:35 PM, sander at dune-project.org wrote:
> Author: sander
> Date: 2010-03-26 22:35:56 +0100 (Fri, 26 Mar 2010)
> New Revision: 6536
>
> Modified:
> trunk/dune/grid/uggrid/uggridintersections.cc
> Log:
> Use the GenericGeometries facilities to compute intersection normals.
>
> This may be a bit slower now, but the code is *a lot* shorter.
> The current implementation sets up a BasicGeometry for the inside element,
> then computes local coordinates with respect to this element and then
> uses the normal-computing method of BasicGeometry. Speed could be
> regained if there was something like a generic IntersectionGeometry,
> which would save me getting corner coordinates I don't need and
> doing the coordinate transformation.
>
>
> Modified: trunk/dune/grid/uggrid/uggridintersections.cc
> ===================================================================
> --- trunk/dune/grid/uggrid/uggridintersections.cc 2010-03-26 21:11:14 UTC (rev 6535)
> +++ trunk/dune/grid/uggrid/uggridintersections.cc 2010-03-26 21:35:56 UTC (rev 6536)
> @@ -5,87 +5,26 @@
>
> #include <set>
>
> -// const Dune::FieldVector<typename GridImp::ctype, Dune::UGGridLevelIntersection::dimworld>& Dune::UGGridLevelIntersection<GridImp>::outerNormal(const Dune::FieldVector<typename GridImp::ctype, (Dune::UGGridLevelIntersection<GridImp>::dim - 1)>&) const;
> -// const Dune::FieldVector<typename GridImp::ctype, Dune::UGGridLevelIntersection<GridImp>::dimworld>& Dune::UGGridLevelIntersection<GridImp>::outerNormal(const Dune::FieldVector<typename GridImp::ctype, (Dune::UGGridLevelIntersection<GridImp>::dim - 1)>&) const;
>
> template<class GridImp>
> const typename Dune::UGGridLevelIntersection<GridImp>::WorldVector&
> Dune::UGGridLevelIntersection<GridImp>::outerNormal
> (const typename Dune::UGGridLevelIntersection<GridImp>::FaceVector& local) const
> {
> - // //////////////////////////////////////////////////////
> - // Implementation for 3D
> - // //////////////////////////////////////////////////////
> + // Implementation uses the GenericGeometry facilities. For this we need to set up a
> + // generic geometry of the _inside_ element, which then has a method computing the normal.
> + std::vector<FieldVector<typename GridImp::ctype, dim> > corners(inside()->template count<dim>());
>
> - if (dim == 3) {
> + for (int i=0; i<corners.size(); i++)
> + corners[i] = inside()->geometry().corner(i);
>
> - if (UG_NS<dim>::Corners_Of_Side(center_, neighborCount_) == 3) {
> + GenericGeometry::BasicGeometry<dim, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,dim,dim> > insideGeometry(inside()->type(), corners);
>
> - // A triangular intersection. The normals are constant
> - const UGCtype* aPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 0))->myvertex->iv.x;
> - const UGCtype* bPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 1))->myvertex->iv.x;
> - const UGCtype* cPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 2))->myvertex->iv.x;
> -
> - FieldVector<UGCtype, 3> ba, ca;
> -
> - for (int i=0; i<3; i++) {
> - ba[i] = bPos[i] - aPos[i];
> - ca[i] = cPos[i] - aPos[i];
> - }
> -
> - outerNormal_[0] = ba[1]*ca[2] - ba[2]*ca[1];
> - outerNormal_[1] = ba[2]*ca[0] - ba[0]*ca[2];
> - outerNormal_[2] = ba[0]*ca[1] - ba[1]*ca[0];
> -
> - } else {
> -
> - // A quadrilateral: compute the normal in each corner and do bilinear interpolation
> - // The cornerNormals array uses UG corner numbering
> - FieldVector<UGCtype,3> cornerNormals[4];
> - for (int i=0; i<4; i++) {
> -
> - // Compute the normal on the i-th corner
> - const UGCtype* aPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_,neighborCount_,i))->myvertex->iv.x;
> - const UGCtype* bPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_,neighborCount_,(i+1)%4))->myvertex->iv.x;
> - const UGCtype* cPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_,neighborCount_,(i+3)%4))->myvertex->iv.x;
> -
> - FieldVector<UGCtype, 3> ba, ca;
> -
> - for (int j=0; j<3; j++) {
> - ba[j] = bPos[j] - aPos[j];
> - ca[j] = cPos[j] - aPos[j];
> - }
> -
> - cornerNormals[i][0] = ba[1]*ca[2] - ba[2]*ca[1];
> - cornerNormals[i][1] = ba[2]*ca[0] - ba[0]*ca[2];
> - cornerNormals[i][2] = ba[0]*ca[1] - ba[1]*ca[0];
> - }
> + // Actually compute the normal.
> + // Note: The local coordinates that have to be provided are with respect to the inside() element
> + outerNormal_ = insideGeometry.normal(indexInInside(),
> + geometryInInside().global(local));
>
> - // Bilinear interpolation
> - for (int i=0; i<3; i++)
> - outerNormal_[i] = (1-local[0])*(1-local[1])*cornerNormals[0][i]
> - + local[0] * (1-local[1]) * cornerNormals[1][i]
> - + local[0] * local[1] * cornerNormals[2][i]
> - + (1-local[0]) * local[1] * cornerNormals[3][i];
> -
> - }
> -
> - } else { // if (dim==3) ... else
> -
> - // //////////////////////////////////////////////////////
> - // Implementation for 2D
> - // //////////////////////////////////////////////////////
> -
> - // Get the vertices of this side.
> - const UGCtype* aPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 0))->myvertex->iv.x;
> - const UGCtype* bPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 1))->myvertex->iv.x;
> -
> - // compute normal
> - outerNormal_[0] = bPos[1] - aPos[1];
> - outerNormal_[1] = aPos[0] - bPos[0];
> -
> - }
> -
> return outerNormal_;
> }
>
> @@ -208,85 +147,30 @@
> // Implementations for the class UGGridLeafIntersection
> // /////////////////////////////////////////////////////////////////////////////
>
> -/** \bug Doesn't work properly for nonplanar nonconforming quadrilateral faces */
> +/** \bug Doesn't work properly for nonplanar nonconforming quadrilateral faces,
> + because the local coordinates are interpreted as being with respect to the element
> + face. Instead, they should be interpreted with respect to the intersection.
> + If the face is flat this doesn't matter.
> +*/
> template<class GridImp>
> const typename Dune::UGGridLeafIntersection<GridImp>::WorldVector&
> Dune::UGGridLeafIntersection<GridImp>::outerNormal
> (const typename Dune::UGGridLeafIntersection<GridImp>::FaceVector& local) const
> {
> - // //////////////////////////////////////////////////////
> - // Implementation for 3D
> - // //////////////////////////////////////////////////////
> + // Implementation uses the GenericGeometry facilities. For this we need to set up a
> + // generic geometry of the _inside_ element, which then has a method computing the normal.
> + std::vector<FieldVector<typename GridImp::ctype, dim> > corners(inside()->template count<dim>());
>
> - if (dim == 3) {
> + for (int i=0; i<corners.size(); i++)
> + corners[i] = inside()->geometry().corner(i);
>
> - if (UG_NS<dim>::Corners_Of_Side(center_, neighborCount_) == 3) {
> + GenericGeometry::BasicGeometry<dim, GenericGeometry::DefaultGeometryTraits<typename GridImp::ctype,dim,dim> > insideGeometry(inside()->type(), corners);
>
> - // A triangular intersection. The normals are constant
> - const UGCtype* aPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 0))->myvertex->iv.x;
> - const UGCtype* bPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 1))->myvertex->iv.x;
> - const UGCtype* cPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 2))->myvertex->iv.x;
> -
> - FieldVector<UGCtype, 3> ba, ca;
> -
> - for (int i=0; i<3; i++) {
> - ba[i] = bPos[i] - aPos[i];
> - ca[i] = cPos[i] - aPos[i];
> - }
> -
> - outerNormal_[0] = ba[1]*ca[2] - ba[2]*ca[1];
> - outerNormal_[1] = ba[2]*ca[0] - ba[0]*ca[2];
> - outerNormal_[2] = ba[0]*ca[1] - ba[1]*ca[0];
> -
> - } else {
> -
> - // A quadrilateral: compute the normal in each corner and do bilinear interpolation
> - // The cornerNormals array uses UG corner numbering
> - FieldVector<UGCtype,3> cornerNormals[4];
> - for (int i=0; i<4; i++) {
> -
> - // Compute the normal on the i-th corner
> - const UGCtype* aPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_,neighborCount_,i))->myvertex->iv.x;
> - const UGCtype* bPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_,neighborCount_,(i+1)%4))->myvertex->iv.x;
> - const UGCtype* cPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_,neighborCount_,(i+3)%4))->myvertex->iv.x;
> -
> - FieldVector<UGCtype, 3> ba, ca;
> -
> - for (int j=0; j<3; j++) {
> - ba[j] = bPos[j] - aPos[j];
> - ca[j] = cPos[j] - aPos[j];
> - }
> -
> - cornerNormals[i][0] = ba[1]*ca[2] - ba[2]*ca[1];
> - cornerNormals[i][1] = ba[2]*ca[0] - ba[0]*ca[2];
> - cornerNormals[i][2] = ba[0]*ca[1] - ba[1]*ca[0];
> - }
> + // Actually compute the normal.
> + // Note: The local coordinates that have to be provided are with respect to the inside() element
> + outerNormal_ = insideGeometry.normal(indexInInside(),
> + geometryInInside().global(local));
>
> - // Bilinear interpolation
> - for (int i=0; i<3; i++)
> - outerNormal_[i] = (1-local[0])*(1-local[1])*cornerNormals[0][i]
> - + local[0] * (1-local[1]) * cornerNormals[1][i]
> - + local[0] * local[1] * cornerNormals[2][i]
> - + (1-local[0]) * local[1] * cornerNormals[3][i];
> -
> - }
> -
> - } else { // if (dim==3) ... else
> -
> - // //////////////////////////////////////////////////////
> - // Implementation for 2D
> - // //////////////////////////////////////////////////////
> -
> - // Get the vertices of this side.
> - const UGCtype* aPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 0))->myvertex->iv.x;
> - const UGCtype* bPos = UG_NS<dim>::Corner(center_,UG_NS<dim>::Corner_Of_Side(center_, neighborCount_, 1))->myvertex->iv.x;
> -
> - // compute normal
> - outerNormal_[0] = bPos[1] - aPos[1];
> - outerNormal_[1] = aPos[0] - bPos[0];
> -
> - }
> -
> return outerNormal_;
> }
>
>
>
> _______________________________________________
> Dune-Commit mailing list
> Dune-Commit at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-commit
--
Martin Nolte <nolte at mathematik.uni-freiburg.de>
Universität Freiburg phone:
+49-761-203-5642
Abteilung für angewandte Mathematik fax:
+49-761-203-5632
Hermann-Herder-Straße 10
79104 Freiburg, Germany
More information about the Dune
mailing list