[Dune] [Dune-Commit] dune-grid r6536 - trunk/dune/grid/uggrid
Oliver Sander
sander at mi.fu-berlin.de
Sat Mar 27 16:42:06 CET 2010
Hi Martin!
Thanks for the information. I know that in general I need the entire
element
to compute its surface normals. My idea was that if I know that
gridDim==worldDim
then I can compute the normal knowing only the face vertices. The more
I think
about it the more this appears to specialized to have it in
GenericGeometries.
Maybe I'll implement something like this in UGGrid directly.
best,
Oliver
Martin Nolte schrieb:
> 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
>>
>
>
--
************************************************************************
* Oliver Sander ** email: sander at mi.fu-berlin.de *
* Freie Universität Berlin ** phone: + 49 (30) 838 75348 *
* Institut für Mathematik ** URL : page.mi.fu-berlin.de/~sander *
* Arnimallee 6 ** -------------------------------------*
* 14195 Berlin, Germany ** Member of MATHEON (www.matheon.de) *
************************************************************************
More information about the Dune
mailing list