[Dune] Fwd: Re: [Dune-Commit] dune-geometry r239 - trunk/dune/geometry/genericgeometry/test

[Martin Nolte]

-------- Original Message --------
Subject: Re: [Dune-Commit] dune-geometry r239 - 
trunk/dune/geometry/genericgeometry/test
Date: Mon, 10 Sep 2012 13:12:13 +0200
From: Martin Nolte <nolte at mathematik.uni-freiburg.de>
To: dune-commit at dune-project.org

Hi Oli,

thanks for the work. I did forget to update the test.

Just one question: What's wrong with

jtAsFieldMatrix = jt;

or

jtAsFieldMatrix
   = static_cast< FieldMatrix< ctype, mydim, coorddim > >( jt );

Both should work for all implementations that I know. For SPGrid, you
currently can even write

const FieldMatrix< ctype, mydim, coorddim > &jtAsFieldMatrix = jt;

Best,

Martin


On 09/10/2012 12:53 PM, sander at dune-project.org wrote:
> Author: sander
> Date: 2012-09-10 12:53:47 +0200 (Mon, 10 Sep 2012)
> New Revision: 239
>
> Modified:
>     trunk/dune/geometry/genericgeometry/test/checkgeometry.hh
> Log:
> Test using the exported types JacobianTransposed and JacobianInverseTransposed
>
> We cannot assume that these are FieldMatrices anymore.
>
> This patch contains some evil trickery.  In the test we do need FieldMatrices,
> e.g., to get operator[] and other niceties.  However I couldn't find a simple
> universal way to get a FieldMatrix from JacobianTransposed.  Therefore the
> code now multiplies the matrix from the right by the identity.  The result
> is a FieldMatrix :-)
>
> Modified: trunk/dune/geometry/genericgeometry/test/checkgeometry.hh
> ===================================================================
> --- trunk/dune/geometry/genericgeometry/test/checkgeometry.hh	2012-09-07 07:53:28 UTC (rev 238)
> +++ trunk/dune/geometry/genericgeometry/test/checkgeometry.hh	2012-09-10 10:53:47 UTC (rev 239)
> @@ -51,6 +51,12 @@
>       // vector type used for image points
>       typedef typename TestGeometry::GlobalCoordinate GlobalCoordinate;
>
> +    // Matrix-like type for the return value of the jacobianTransposed method
> +    typedef typename TestGeometry::JacobianTransposed JacobianTransposed;
> +
> +    // Matrix-like type for the return value of the jacobianInverseTransposed method
> +    typedef typename TestGeometry::JacobianInverseTransposed JacobianInverseTransposed;
> +
>       ////////////////////////////////////////////////////////////////////////
>
>       GeometryType type = geometry.type();
> @@ -101,11 +107,45 @@
>
>         // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed'
>         // return matrices that are inverse to each other.
> -      const FieldMatrix< ctype, mydim, coorddim > &jt = geometry.jacobianTransposed( x );
> -      const FieldMatrix< ctype, coorddim, mydim > &jit = geometry.jacobianInverseTransposed( x );
> +      const JacobianTransposed         &jt = geometry.jacobianTransposed( x );
> +      const JacobianInverseTransposed &jit = geometry.jacobianInverseTransposed( x );
>
> +      // Transform to FieldMatrix, so we can have coefficent access and other goodies
> +      // We need some black magic for the transformation, because there is no
> +      // official easy way yet.
> +      // The following code does the transformation by multiplying jt and jit from
> +      // the right by identity matrices.  That way, only the mv method is used.
> +      FieldMatrix< ctype, mydim, coorddim > jtAsFieldMatrix;
> +      for (int j=0; j<coorddim; j++) {
> +
> +          FieldVector<ctype,coorddim> idColumn(0);
> +          idColumn[j] = 1;
> +
> +          FieldVector<ctype,mydim> column;
> +          jt.mv(idColumn,column);
> +
> +          for (int k=0; k<mydim; k++)
> +              jtAsFieldMatrix[k][j] = column[k];
> +
> +      }
> +
> +      FieldMatrix< ctype, coorddim, mydim > jitAsFieldMatrix;
> +      for (int j=0; j<mydim; j++) {
> +
> +          FieldVector<ctype,mydim> idColumn(0);
> +          idColumn[j] = 1;
> +
> +          FieldVector<ctype,coorddim> column;
> +          jit.mv(idColumn,column);
> +
> +          for (int k=0; k<coorddim; k++)
> +              jitAsFieldMatrix[k][j] = column[k];
> +
> +      }
> +
> +
>         FieldMatrix< ctype, mydim, mydim > id;
> -      FMatrixHelp::multMatrix( jt, jit, id );
> +      FMatrixHelp::multMatrix( jtAsFieldMatrix, jitAsFieldMatrix, id );
>         bool isId = true;
>         for( int j = 0; j < mydim; ++j )
>           for( int k = 0; k < mydim; ++k )
> @@ -130,7 +170,8 @@
>         for( int i = 0; i < mydim; ++i )
>           for( int j = 0; j < mydim; ++j )
>             for( int k = 0; k < coorddim; ++k )
> -            jtj[ i ][ j ] += jt[ i ][ k ] * jt[ j ][ k ];
> +            jtj[ i ][ j ] += jtAsFieldMatrix[ i ][ k ] * jtAsFieldMatrix[ j ][ k ];
> +
>         if( std::abs( std::sqrt( jtj.determinant() ) - geometry.integrationElement( x ) ) > 1e-8 ) {
>           std::cerr << "Error: integrationElement is not consistent with jacobianTransposed." << std::endl;
>           pass = false;
>
>
> _______________________________________________
> Dune-Commit mailing list
> Dune-Commit at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-commit
>

-- 
Dr. Martin Nolte <nolte at mathematik.uni-freiburg.de>

Universität Freiburg                                   phone: +49-761-203-5630
Abteilung für angewandte Mathematik                    fax:   +49-761-203-5632
Hermann-Herder-Straße 10
79104 Freiburg, Germany