← Back to branch summary

~ubuntu-branches/ubuntu/utopic/dune-grid/utopic-proposed

« back to all changes in this revision

Viewing changes to dune/grid/alugrid/2d/geometry_imp.cc

  • Committer: Package Import Robot
  • Author(s): Ansgar Burchardt
  • Date: 2014-02-14 10:49:16 UTC
  • mfrom: (1.3.1) (5.1.4 experimental)
  • Revision ID: package-import@ubuntu.com-20140214104916-2clchnny3eu7ks4w
Tags: 2.3.0-2
Install /usr/share/dune-grid.

Show diffs side-by-side

added added

removed removed

Lines of Context:
dune/grid/alugrid/2d/geometry_imp.cc
 
1
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 
2
// vi: set et ts=4 sw=2 sts=2:
1
3
#ifndef DUNE_ALU2DGRID_GEOMETRYIMP_CC
2
4
#define DUNE_ALU2DGRID_GEOMETRYIMP_CC
3
5
 
4
 
#include <dune/geometry/genericgeometry/conversion.hh>
5
6
#include <dune/geometry/referenceelements.hh>
6
7
 
7
8
namespace Dune
8
9
{
9
10
 
10
 
// implementation of methods
11
 
 
12
 
//**********************************************************************
13
 
//
14
 
// --ALU2dGridGeometry
15
 
// --Geometry
16
 
//**********************************************************************
17
 
 
18
 
 
19
 
template< int mydim, int cdim, class GridImp >
20
 
inline ALU2dGridGeometry< mydim, cdim, GridImp >::ALU2dGridGeometry ()
21
 
: geoImpl_(),
22
 
  det_( 1.0 )
23
 
{}
24
 
 
25
 
 
26
 
//! print the GeometryInformation
27
 
template< int mydim, int cdim, class GridImp >
28
 
inline void ALU2dGridGeometry< mydim, cdim, GridImp >::print ( std::ostream &out ) const
29
 
{
30
 
  out << "ALU2dGridGeometry< " << mydim << ", " << cdim << " > = ";
31
 
  char c = '{';
32
 
  for( int i = 0; i < corners(); ++i )
33
 
  {
34
 
    out << c << " (" << corner( i ) << ")";
35
 
    c = ',';
36
 
  }
37
 
  out << " }" << std::endl;
38
 
}
39
 
 
40
 
 
41
 
///////////////////////////////////////////////////////////////////////
42
 
 
43
 
 
44
 
//! access to coordinates of corners. Index is the number of the corner 
45
 
template <int mydim, int cdim, class GridImp>
46
 
inline const typename ALU2dGridGeometry< mydim, cdim, GridImp >::GlobalCoordinate &
47
 
ALU2dGridGeometry< mydim, cdim, GridImp >::operator[] ( int i ) const
48
 
{
49
 
  typedef GenericGeometry::MapNumberingProvider< mydim > Numbering;
50
 
  const unsigned int tid = type().id();
51
 
  static GlobalCoordinate c;
52
 
  return (c = corner( Numbering::template dune2generic< mydim >( tid, i ) ));
53
 
}
54
 
 
55
 
 
56
 
//! access to coordinates of corners. Index is the number of the corner 
57
 
template< int mydim, int cdim, class GridImp >
58
 
inline typename ALU2dGridGeometry< mydim, cdim, GridImp >::GlobalCoordinate
59
 
ALU2dGridGeometry< mydim, cdim, GridImp >::corner ( int i ) const
60
 
{
61
 
  const GenericReferenceElement< alu2d_ctype, mydim > &refElement
62
 
    = GenericReferenceElements< alu2d_ctype, mydim >::general( type() );
63
 
  return global( refElement.position( i, mydim ) );
64
 
}
65
 
 
66
 
 
67
 
//! maps a local coordinate within reference element to 
68
 
//! global coordinate in element 
69
 
template< int mydim, int cdim, class GridImp >
70
 
inline typename ALU2dGridGeometry< mydim, cdim, GridImp >::GlobalCoordinate
71
 
ALU2dGridGeometry< mydim, cdim, GridImp >::global ( const LocalCoordinate &local ) const
72
 
{
73
 
  GlobalCoordinate global;
74
 
  geoImpl_.map2world( local, global );
75
 
  return global;
76
 
}
77
 
 
78
 
 
79
 
//! maps a global coordinate within the element to a 
80
 
//! local coordinate in its reference element
81
 
template< int mydim, int cdim, class GridImp >
82
 
inline typename ALU2dGridGeometry< mydim, cdim, GridImp >::LocalCoordinate
83
 
ALU2dGridGeometry< mydim, cdim, GridImp >::local ( const GlobalCoordinate &global ) const
84
 
{
85
 
  if( mydim == 0 )
86
 
    return LocalCoordinate( 1 );
87
 
 
88
 
  LocalCoordinate local;
89
 
  geoImpl_.world2map( global, local );
90
 
  return local;
91
 
}
92
 
 
93
 
 
94
 
template< int mydim, int cdim, class GridImp >
95
 
inline alu2d_ctype
96
 
ALU2dGridGeometry< mydim, cdim, GridImp >::integrationElement ( const LocalCoordinate &local ) const
97
 
{
98
 
  if ( eltype == ALU2DSPACE triangle || mydim < 2 )
99
 
  {
100
 
    assert( geoImpl_.valid() );
101
 
    return (mydim == 0 ? 1.0 : det_);
102
 
  }
103
 
  else
104
 
    return geoImpl_.det(local);
105
 
}
106
 
 
107
 
 
108
 
template< int mydim, int cdim, class GridImp >
109
 
inline alu2d_ctype ALU2dGridGeometry< mydim, cdim, GridImp >::volume () const
110
 
{
111
 
  assert( geoImpl_.valid() );
112
 
  if( mydim == 2 )
113
 
  {
114
 
    switch( GridImp::elementType )
115
 
    {
116
 
    case ALU2DSPACE triangle:
117
 
      return 0.5 * det_;
118
 
 
119
 
    case ALU2DSPACE quadrilateral:
120
 
      return det_;
121
 
 
122
 
    case ALU2DSPACE mixed:
123
 
      DUNE_THROW( NotImplemented, "Geometry::volume() not implemented for ElementType mixed." );
124
 
    }
125
 
  }
126
 
  else
127
 
    return (mydim == 0 ? 1.0 : det_);
128
 
}
129
 
 
130
 
 
131
 
template< int mydim, int cdim, class GridImp >
132
 
inline const FieldMatrix<alu2d_ctype,mydim,cdim>&
133
 
ALU2dGridGeometry< mydim, cdim, GridImp>::jacobianTransposed ( const LocalCoordinate &local ) const
134
 
{
135
 
  return geoImpl_.jacobianTransposed( local );
136
 
}
137
 
 
138
 
 
139
 
template< int mydim, int cdim, class GridImp >
140
 
inline const FieldMatrix<alu2d_ctype,cdim,mydim>&
141
 
ALU2dGridGeometry< mydim, cdim, GridImp >::jacobianInverseTransposed ( const LocalCoordinate &local ) const
142
 
{
143
 
  return geoImpl_.jacobianInverseTransposed( local );
144
 
}
145
 
 
146
 
 
147
 
//! built Geometry for triangles 
148
 
template< int mydim, int cdim, class GridImp >
149
 
inline bool
150
 
ALU2dGridGeometry< mydim, cdim, GridImp >::buildGeom( const HElementType &item )
151
 
{
152
 
  // check item 
153
 
  assert( &item );
154
 
  assert( mydim == 2 );
155
 
 
156
 
  // update geometry impl 
157
 
  geoImpl_.update( item );
158
 
 
159
 
  // store volume
160
 
  det_ = (geoImpl_.corners() == 3 ? 2 * item.area() : item.area());
161
 
 
162
 
  // geometry built
163
 
  return true;
164
 
}
165
 
 
166
 
 
167
 
//! built Geometry for edges 
168
 
template <int mydim, int cdim, class GridImp>
169
 
inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
170
 
buildGeom(const HElementType & item, const int aluFace)
171
 
{
172
 
  assert( mydim == 1 );
173
 
  // face 1 is twisted 
174
 
  const int nf = item.numfaces();
175
 
 
176
 
  // for triangles face 1 is twisted and for quatrilaterals faces 1 and 2
177
 
  const int twist = (nf == 3) ? (aluFace % 2) : (aluFace>>1)^(aluFace&1);
178
 
  // std::cout << nf << " " << aluFace << " " << twist << std::endl;
179
 
  //           << " ( " item.getVertex( (aluFace + 1 + twist ) % nf )->coord() ,
180
 
  //           << " , " item.getVertex( (aluFace + 2 - twist ) % nf )->coord() 
181
 
  //           << " ) " << std::endl;
182
 
 
183
 
  // check item 
184
 
  assert( &item );
185
 
 
186
 
 
187
 
  // update geometry impl 
188
 
  geoImpl_.update( item.getVertex( (aluFace + 1 + twist ) % nf )->coord() ,
189
 
                   item.getVertex( (aluFace + 2 - twist ) % nf )->coord() );
190
 
 
191
 
  // store volume
192
 
  det_ = item.sidelength( aluFace );
193
 
  //assert( std::abs( det_ - geoImpl_.det( LocalCoordinate(0.5) ) ) < 1e-14 );
194
 
  
195
 
  // geometry built
196
 
  return true;
197
 
}
198
 
 
199
 
//! built Geometry for vertices 
200
 
template <int mydim, int cdim, class GridImp>
201
 
inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
202
 
buildGeom(const VertexType & item , const int )
203
 
{
204
 
  assert( mydim == 0 );
205
 
 
206
 
  assert( &item );
207
 
  // update geometry impl 
208
 
  geoImpl_.update( item.coord() );
209
 
 
210
 
  // volume is already 1.0 
211
 
 
212
 
  return true; 
213
 
}
214
 
 
215
 
// built Geometry
216
 
template <int mydim, int cdim, class GridImp>
217
 
template <class GeometryType, class LocalGeometryType >
218
 
inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
219
 
buildLocalGeom(const GeometryType &geo, const LocalGeometryType & localGeom)
220
 
{
221
 
  // update geometry 
222
 
  geoImpl_.updateLocal( geo, localGeom );
223
 
 
224
 
  // calculate volume
225
 
  LocalCoordinate local( 0.25 );
226
 
  det_ = geoImpl_.det( local );
227
 
  
228
 
  // geometry built
229
 
  return true;
230
 
}
231
 
 
232
 
// built Geometry (faceNumber is in generic numbering)
233
 
template< int mydim, int cdim, class GridImp >
234
 
inline std::pair< FieldMatrix< alu2d_ctype, 4, 2 >, FieldVector< alu2d_ctype, 4 > >
235
 
ALU2dGridGeometry< mydim, cdim, GridImp >::calculateReferenceCoords ( const int corners )
236
 
{
237
 
  // calculate reference coordinates of aLUGrid reference triangle
238
 
 
239
 
  FieldMatrix< alu2d_ctype, 4, 2 > refCoord( 0. );
240
 
  FieldVector< alu2d_ctype, 4 > lengths( 1. );
241
 
 
242
 
  // point 1 
243
 
  refCoord[1][0] = 1.0;
244
 
  // point (corners-1)
245
 
  refCoord[corners-1][1] = 1.0;
246
 
  if( corners == 3 )
247
 
  {
248
 
    lengths[0] = M_SQRT2;
249
 
  } 
250
 
  else
251
 
  {
252
 
    // point 2
253
 
    refCoord[2][0] = 1.0;
254
 
    refCoord[2][1] = 1.0;
255
 
  }
256
 
  return std::make_pair( refCoord, lengths );
257
 
}
258
 
 
259
 
// built Geometry (faceNumber is in generic numbering)
260
 
template <int mydim, int cdim, class GridImp>
261
 
inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
262
 
buildLocalGeometry(const int aluFace, const int twist, const int corners)
263
 
{
264
 
  assert( twist == 0 || twist == 1 );
265
 
  assert( mydim == 1 );
266
 
 
267
 
  // get coordinates of reference element
268
 
  typedef std::pair< FieldMatrix< alu2d_ctype, 4, 2 >, FieldVector< alu2d_ctype, 4 > > RefCoord;
269
 
  RefCoord refCoord( calculateReferenceCoords( corners ) );
270
 
 
271
 
  geoImpl_.update( refCoord.first[ ( aluFace + 1+twist ) % corners ],
272
 
                   refCoord.first[ ( aluFace + 2-twist ) % corners ] );
273
 
 
274
 
  // get length of faces 
275
 
  det_ = refCoord.second[ aluFace ];
276
 
  
277
 
  // geometry built
278
 
  return true;
279
 
}
280
 
 
281
 
// built Geometry
282
 
template <int mydim, int cdim, class GridImp >
283
 
inline bool ALU2dGridGeometry<mydim, cdim,GridImp>::
284
 
buildGeomInFather(const Geometry & fatherGeom , 
285
 
                  const Geometry & myGeom)
286
 
{
287
 
  // update geometry 
288
 
  geoImpl_.updateLocal( fatherGeom, myGeom );
289
 
 
290
 
  // store volume which is a part of one 
291
 
  det_ = myGeom.volume() / fatherGeom.volume();
292
 
  assert( (det_ > 0.0) && (det_ < 1.0) );
293
 
 
294
 
  return true;
295
 
}
296
 
  
297
 
} //end namespace Dune 
 
11
  // implementation of methods
 
12
 
 
13
  //**********************************************************************
 
14
  //
 
15
  // --ALU2dGridGeometry
 
16
  // --Geometry
 
17
  //**********************************************************************
 
18
 
 
19
 
 
20
  template< int mydim, int cdim, class GridImp >
 
21
  inline ALU2dGridGeometry< mydim, cdim, GridImp >::ALU2dGridGeometry ()
 
22
  {
 
23
    getObject();
 
24
  }
 
25
 
 
26
  template< int mydim, int cdim, class GridImp >
 
27
  inline ALU2dGridGeometry< mydim, cdim, GridImp >&
 
28
  ALU2dGridGeometry< mydim, cdim, GridImp >::operator = ( const ALU2dGridGeometry& other )
 
29
  {
 
30
    if( &other != this )
 
31
    {
 
32
      removeObj();
 
33
      assign( other );
 
34
    }
 
35
    return *this;
 
36
  }
 
37
 
 
38
  template< int mydim, int cdim, class GridImp >
 
39
  inline ALU2dGridGeometry< mydim, cdim, GridImp >::ALU2dGridGeometry ( const ALU2dGridGeometry& other )
 
40
  {
 
41
    assign( other );
 
42
  }
 
43
 
 
44
  template< int mydim, int cdim, class GridImp >
 
45
  inline ALU2dGridGeometry< mydim, cdim, GridImp >::~ALU2dGridGeometry ()
 
46
  {
 
47
    removeObj();
 
48
  }
 
49
 
 
50
  template< int mydim, int cdim, class GridImp >
 
51
  inline void ALU2dGridGeometry< mydim, cdim, GridImp >::getObject()
 
52
  {
 
53
    geoImpl_ = geoProvider().getEmptyObject();
 
54
    geoImpl().reset();
 
55
  }
 
56
 
 
57
  template< int mydim, int cdim, class GridImp >
 
58
  inline void ALU2dGridGeometry< mydim, cdim, GridImp >::assign( const ALU2dGridGeometry& other )
 
59
  {
 
60
    // copy pointer
 
61
    geoImpl_ = other.geoImpl_ ;
 
62
 
 
63
    // increase reference count
 
64
    ++ geoImpl();
 
65
  }
 
66
 
 
67
  template< int mydim, int cdim, class GridImp >
 
68
  inline void ALU2dGridGeometry< mydim, cdim, GridImp >::removeObj()
 
69
  {
 
70
    // decrease reference count
 
71
    -- geoImpl();
 
72
 
 
73
    // if reference count is zero free the object
 
74
    if( ! geoImpl() )
 
75
    {
 
76
      geoProvider().freeObject( geoImpl_ );
 
77
    }
 
78
 
 
79
    // reset pointer
 
80
    geoImpl_ = 0;
 
81
  }
 
82
 
 
83
  template< int mydim, int cdim, class GridImp >
 
84
  inline void ALU2dGridGeometry< mydim, cdim, GridImp >::invalidate()
 
85
  {
 
86
    // if geometry is used elsewhere remove the pointer
 
87
    // and get a new one
 
88
    if( geoImpl().stillUsed() )
 
89
    {
 
90
      // remove old object
 
91
      removeObj();
 
92
 
 
93
      // get new object
 
94
      getObject();
 
95
    }
 
96
    else
 
97
    {
 
98
      // otherwise invalidate object
 
99
      geoImpl().invalidate();
 
100
    }
 
101
  }
 
102
 
 
103
  //! print the GeometryInformation
 
104
  template< int mydim, int cdim, class GridImp >
 
105
  inline void ALU2dGridGeometry< mydim, cdim, GridImp >::print ( std::ostream &out ) const
 
106
  {
 
107
    out << "ALU2dGridGeometry< " << mydim << ", " << cdim << " > = ";
 
108
    char c = '{';
 
109
    for( int i = 0; i < corners(); ++i )
 
110
    {
 
111
      out << c << " (" << corner( i ) << ")";
 
112
      c = ',';
 
113
    }
 
114
    out << " }" << std::endl;
 
115
  }
 
116
 
 
117
 
 
118
  ///////////////////////////////////////////////////////////////////////
 
119
 
 
120
 
 
121
  //! access to coordinates of corners. Index is the number of the corner
 
122
  template< int mydim, int cdim, class GridImp >
 
123
  inline typename ALU2dGridGeometry< mydim, cdim, GridImp >::GlobalCoordinate
 
124
  ALU2dGridGeometry< mydim, cdim, GridImp >::corner ( int i ) const
 
125
  {
 
126
    return geoImpl().corner( i );
 
127
  }
 
128
 
 
129
 
 
130
  //! maps a local coordinate within reference element to
 
131
  //! global coordinate in element
 
132
  template< int mydim, int cdim, class GridImp >
 
133
  inline typename ALU2dGridGeometry< mydim, cdim, GridImp >::GlobalCoordinate
 
134
  ALU2dGridGeometry< mydim, cdim, GridImp >::global ( const LocalCoordinate &local ) const
 
135
  {
 
136
    GlobalCoordinate global;
 
137
    geoImpl().map2world( local, global );
 
138
    return global;
 
139
  }
 
140
 
 
141
 
 
142
  //! maps a global coordinate within the element to a
 
143
  //! local coordinate in its reference element
 
144
  template< int mydim, int cdim, class GridImp >
 
145
  inline typename ALU2dGridGeometry< mydim, cdim, GridImp >::LocalCoordinate
 
146
  ALU2dGridGeometry< mydim, cdim, GridImp >::local ( const GlobalCoordinate &global ) const
 
147
  {
 
148
    if( mydim == 0 )
 
149
      return LocalCoordinate( 1 );
 
150
 
 
151
    LocalCoordinate local;
 
152
    geoImpl().world2map( global, local );
 
153
    return local;
 
154
  }
 
155
 
 
156
 
 
157
  template< int mydim, int cdim, class GridImp >
 
158
  inline alu2d_ctype
 
159
  ALU2dGridGeometry< mydim, cdim, GridImp >::integrationElement ( const LocalCoordinate &local ) const
 
160
  {
 
161
    assert( geoImpl().valid() );
 
162
    if( mydim == cdim )
 
163
    {
 
164
      if( std::abs( geoImpl().det( local ) - jacobianTransposed( local
 
165
                                                                 ).determinant() ) > 1e-8 )
 
166
        std::cout << "det = " << geoImpl().det( local ) << "  inv.det " <<
 
167
        jacobianTransposed( local ).determinant() << std::endl;
 
168
      assert( std::abs( geoImpl().det( local ) - jacobianTransposed( local
 
169
                                                                     ).determinant() ) < 1e-8 );
 
170
    }
 
171
    return geoImpl().det( local );
 
172
  }
 
173
 
 
174
 
 
175
  template< int mydim, int cdim, class GridImp >
 
176
  inline alu2d_ctype ALU2dGridGeometry< mydim, cdim, GridImp >::volume () const
 
177
  {
 
178
    assert( geoImpl().valid() );
 
179
    return geoImpl().volume();
 
180
  }
 
181
 
 
182
 
 
183
  template< int mydim, int cdim, class GridImp >
 
184
  inline const FieldMatrix<alu2d_ctype,mydim,cdim>&
 
185
  ALU2dGridGeometry< mydim, cdim, GridImp>::jacobianTransposed ( const LocalCoordinate &local ) const
 
186
  {
 
187
    return geoImpl().jacobianTransposed( local );
 
188
  }
 
189
 
 
190
 
 
191
  template< int mydim, int cdim, class GridImp >
 
192
  inline const FieldMatrix<alu2d_ctype,cdim,mydim>&
 
193
  ALU2dGridGeometry< mydim, cdim, GridImp >::jacobianInverseTransposed ( const LocalCoordinate &local ) const
 
194
  {
 
195
    return geoImpl().jacobianInverseTransposed( local );
 
196
  }
 
197
 
 
198
 
 
199
  //! built Geometry for triangles
 
200
  template< int mydim, int cdim, class GridImp >
 
201
  inline bool
 
202
  ALU2dGridGeometry< mydim, cdim, GridImp >::buildGeom( const HElementType &item )
 
203
  {
 
204
    // check item
 
205
    assert( &item );
 
206
    assert( mydim == 2 );
 
207
 
 
208
    // update geometry impl
 
209
    geoImpl().update( item );
 
210
 
 
211
    // geometry built
 
212
    return true;
 
213
  }
 
214
 
 
215
 
 
216
  //! built Geometry for edges
 
217
  template <int mydim, int cdim, class GridImp>
 
218
  inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
 
219
  buildGeom(const HElementType & item, const int aluFace)
 
220
  {
 
221
    assert( mydim == 1 );
 
222
    // face 1 is twisted
 
223
    const int nf = item.numfaces();
 
224
 
 
225
    // for triangles face 1 is twisted and for quatrilaterals faces 1 and 2
 
226
    const int twist = (nf == 3) ? (aluFace % 2) : (aluFace>>1)^(aluFace&1);
 
227
    // std::cout << nf << " " << aluFace << " " << twist << std::endl;
 
228
    //           << " ( " item.getVertex( (aluFace + 1 + twist ) % nf )->coord() ,
 
229
    //           << " , " item.getVertex( (aluFace + 2 - twist ) % nf )->coord()
 
230
    //           << " ) " << std::endl;
 
231
 
 
232
    // check item
 
233
    assert( &item );
 
234
 
 
235
 
 
236
    // update geometry impl
 
237
    geoImpl().update( item.getVertex( (aluFace + 1 + twist ) % nf )->coord() ,
 
238
                      item.getVertex( (aluFace + 2 - twist ) % nf )->coord() ,
 
239
                      item.sidelength( aluFace )
 
240
                      );
 
241
 
 
242
    // geometry built
 
243
    return true;
 
244
  }
 
245
 
 
246
  //! built Geometry for vertices
 
247
  template <int mydim, int cdim, class GridImp>
 
248
  inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
 
249
  buildGeom(const VertexType & item , const int )
 
250
  {
 
251
    assert( mydim == 0 );
 
252
 
 
253
    assert( &item );
 
254
    // update geometry impl
 
255
    // volume is already 1.0
 
256
    geoImpl().update( item.coord() );
 
257
 
 
258
    // geometry built
 
259
    return true;
 
260
  }
 
261
 
 
262
  // built Geometry
 
263
  template <int mydim, int cdim, class GridImp>
 
264
  template <class GeometryType, class LocalGeometryType >
 
265
  inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
 
266
  buildLocalGeom(const GeometryType &geo, const LocalGeometryType & localGeom)
 
267
  {
 
268
    // update geometry
 
269
    geoImpl().updateLocal( geo, localGeom );
 
270
    // geometry built
 
271
    return true;
 
272
  }
 
273
 
 
274
  // built Geometry (faceNumber is in generic numbering)
 
275
  template< int mydim, int cdim, class GridImp >
 
276
  inline std::pair< FieldMatrix< alu2d_ctype, 4, 2 >, FieldVector< alu2d_ctype, 4 > >
 
277
  ALU2dGridGeometry< mydim, cdim, GridImp >::calculateReferenceCoords ( const int corners )
 
278
  {
 
279
    // calculate reference coordinates of aLUGrid reference triangle
 
280
 
 
281
    FieldMatrix< alu2d_ctype, 4, 2 > refCoord( 0. );
 
282
    FieldVector< alu2d_ctype, 4 > lengths( 1. );
 
283
 
 
284
    // point 1
 
285
    refCoord[1][0] = 1.0;
 
286
    // point (corners-1)
 
287
    refCoord[corners-1][1] = 1.0;
 
288
    if( corners == 3 )
 
289
    {
 
290
      lengths[0] = M_SQRT2;
 
291
    }
 
292
    else
 
293
    {
 
294
      // point 2
 
295
      refCoord[2][0] = 1.0;
 
296
      refCoord[2][1] = 1.0;
 
297
    }
 
298
    return std::make_pair( refCoord, lengths );
 
299
  }
 
300
 
 
301
  // built Geometry (faceNumber is in generic numbering)
 
302
  template <int mydim, int cdim, class GridImp>
 
303
  inline bool ALU2dGridGeometry<mydim,cdim,GridImp>::
 
304
  buildLocalGeometry(const int aluFace, const int twist, const int corners)
 
305
  {
 
306
    assert( twist == 0 || twist == 1 );
 
307
    assert( mydim == 1 );
 
308
 
 
309
    // get coordinates of reference element
 
310
    typedef std::pair< FieldMatrix< alu2d_ctype, 4, 2 >, FieldVector< alu2d_ctype, 4 > > RefCoord;
 
311
    RefCoord refCoord( calculateReferenceCoords( corners ) );
 
312
 
 
313
    geoImpl().update( refCoord.first[ ( aluFace + 1+twist ) % corners ],
 
314
                      refCoord.first[ ( aluFace + 2-twist ) % corners ],
 
315
                      refCoord.second[ aluFace ]
 
316
                      );
 
317
 
 
318
    // geometry built
 
319
    return true;
 
320
  }
 
321
 
 
322
  // built Geometry
 
323
  template <int mydim, int cdim, class GridImp >
 
324
  inline bool ALU2dGridGeometry<mydim, cdim,GridImp>::
 
325
  buildGeomInFather(const Geometry & fatherGeom ,
 
326
                    const Geometry & myGeom)
 
327
  {
 
328
    // update geometry
 
329
    geoImpl().updateLocal( fatherGeom, myGeom );
 
330
    // geometry built
 
331
    return true;
 
332
  }
 
333
 
 
334
} //end namespace Dune
298
335
 
299
336
#endif // #ifndef DUNE_ALU2DGRID_GEOMETRYIMP_CC