SurfaceOperator.h 6.33 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file SurfaceOperator.h */

#ifndef AMDIS_SURFACEOPERATOR_H
#define AMDIS_SURFACEOPERATOR_H

#include "FiniteElemSpace.h"
#include "Mesh.h"
#include "Operator.h"
#include "SurfaceQuadrature.h"
29
#include "OpenMP.h"
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

namespace AMDiS {

  /** 
   * \ingroup Integration
   *
   * \brief
   * A SurfaceOperator is a Operator used for surface integration. Instead of a 
   * normal Quadrature it uses a SurfaceQuadrature with integration points at
   * one element side instead of integration points in the inner of the element.
   * The determinant in ElInfo is replaced by the surface determinant, so the
   * assemblage can be done with the standard Assembler classes.
   * The SurfaceQuadrature is used for the implementation of Neumann and Robin
   * boundary conditions.
   */
  class SurfaceOperator : public Operator
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
48
    /// Creates a SurfaceOperator conforming to operat for the given \ref coords.
49
50
51
52
53
54
55
56
57
    SurfaceOperator(Operator *operat, 
		    VectorOfFixVecs<DimVec<double> > &coords) 
      : Operator(*operat), 
	coords_(coords),
	quad2(NULL),
	quad1GrdPsi(NULL),
	quad1GrdPhi(NULL),
	quad0(NULL)
    {
58
      assembler[omp_get_thread_num()] = NULL;
59
60
61

      int dim = rowFESpace->getMesh()->getDim();
      int degree;
62
      int myRank = omp_get_thread_num();
63
64

      // create surface quadratures
65
      if (secondOrder[myRank].size() > 0) {
66
67
	degree = getQuadratureDegree(2);
	// quad2 = SurfaceQuadrature::provideSurfaceQuadrature(dim, degree, coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
68
	quad2 = new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree),
69
70
71
				      coords);
      }

72
      if (firstOrderGrdPsi[myRank].size() > 0) {
73
	degree = getQuadratureDegree(1, GRD_PSI);
Thomas Witkowski's avatar
Thomas Witkowski committed
74
	quad1GrdPsi = new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree),
75
76
77
78
					    coords);
	//SurfaceQuadrature::provideSurfaceQuadrature(dim, degree, coords);
      }

79
      if (firstOrderGrdPhi[myRank].size() > 0) {
80
	degree = getQuadratureDegree(1, GRD_PHI);
Thomas Witkowski's avatar
Thomas Witkowski committed
81
	quad1GrdPhi = new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree),
82
83
84
85
					    coords);
	//SurfaceQuadrature::provideSurfaceQuadrature(dim, degree, coords);
      }

86
      if (zeroOrder[myRank].size() > 0) {
87
	degree = getQuadratureDegree(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
88
	quad0 = new SurfaceQuadrature(Quadrature::provideQuadrature(dim - 1, degree),
89
90
91
92
93
94
				      coords);
	//SurfaceQuadrature::provideSurfaceQuadrature(dim, degree, coords);
      }

      // initialize assembler with surface quadratures
      optimized = false;
95
      initAssembler(myRank, quad2, quad1GrdPsi, quad1GrdPhi, quad0);
Thomas Witkowski's avatar
Thomas Witkowski committed
96
    }
97

Thomas Witkowski's avatar
Thomas Witkowski committed
98
    /// Adapt surface quadratures to \ref coords.
99
100
101
102
    void adaptSurfaceOperator(VectorOfFixVecs<DimVec<double> > &coords)
    {
      coords_ = coords;

103
      if (quad2)
104
	quad2->scaleSurfaceQuadrature(coords);
105
      if (quad1GrdPsi)
106
	quad1GrdPsi->scaleSurfaceQuadrature(coords);
107
      if (quad1GrdPhi)
108
	quad1GrdPhi->scaleSurfaceQuadrature(coords);
109
      if (quad0)
110
	quad0->scaleSurfaceQuadrature(coords);
Thomas Witkowski's avatar
Thomas Witkowski committed
111
    }
112
113
114
115
116
117
118

    /** \brief
     * Implementation of \ref Operator::getElementMatrix(). Repalces the 
     * determinant by the surface determinant and deligates the call to
     * the base class function.
     */
    virtual void getElementMatrix(const ElInfo *elInfo, 
119
				  ElementMatrix& userMat, 
120
121
122
123
124
125
126
127
				  double factor = 1.0)
    {
      int dim = rowFESpace->getMesh()->getDim();
      double origDet = elInfo->getDet();

      FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT);

      // transform barycentric coords to world coords
128
      for (int i = 0; i < dim; i++)
129
	elInfo->coordToWorld(coords_[i], worldCoords[i]);
130
131

      // set determinant for world coords of the side
Thomas Witkowski's avatar
Thomas Witkowski committed
132
      const_cast<ElInfo*>(elInfo)->setDet(elInfo->calcDet(worldCoords));
133
134
135
136
137
138

      // calc element matrix
      Operator::getElementMatrix(elInfo, userMat, factor);

      // set determinant for world coords of the side
      const_cast<ElInfo*>(elInfo)->setDet(origDet);
Thomas Witkowski's avatar
Thomas Witkowski committed
139
    }
140
141
142
143
144
145
146
  
    /** \brief
     * Implementation of \ref Operator::getElementVector(). Repalces the 
     * determinant by the surface determinant and deligates the call to
     * the base class function.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
147
				  ElementVector& userVec, 
148
149
150
151
152
153
154
155
				  double factor = 1.0) 
    {
      int dim = rowFESpace->getMesh()->getDim();
      double origDet = elInfo->getDet();

      FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT);

      // transform barycentric coords to world coords
156
      for (int i = 0; i < dim; i++)
157
	elInfo->coordToWorld(coords_[i], worldCoords[i]);
158
159

      // set determinant for world coords of the side
Thomas Witkowski's avatar
Thomas Witkowski committed
160
      const_cast<ElInfo*>(elInfo)->setDet(elInfo->calcDet(worldCoords));
161
162
163
164
165

      // calc element vector
      Operator::getElementVector(elInfo, userVec, factor);

      const_cast<ElInfo*>(elInfo)->setDet(origDet);
Thomas Witkowski's avatar
Thomas Witkowski committed
166
    }
167
168
169
170

  protected:
    VectorOfFixVecs<DimVec<double> > coords_;

171
    /// Surface quadratures
172
173
174
175
176
177
178
179
180
    SurfaceQuadrature *quad2;
    SurfaceQuadrature *quad1GrdPsi;
    SurfaceQuadrature *quad1GrdPhi;
    SurfaceQuadrature *quad0;
  };

}

#endif