PenaltyOperator.h 3.55 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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
/** \file PenaltyOperator.h */

#ifndef AMDIS_PENALTYOPERATOR_H
#define AMDIS_PENALTYOPERATOR_H

#include "ElementLevelSet.h"
#include "Flag.h"
#include "Operator.h"
#include "SurfaceOperator.h"
#include "ElementLevelSet.h"

namespace AMDiS {
    class ElInfo;
    class FiniteElemSpace;
}

using namespace AMDiS;
using namespace std;

// ===========================================================================
// ===== class PenaltyOperator ===============================================
// ===========================================================================
//
// Class Description:
//
// ===========================================================================
class PenaltyOperator : public Operator
{
public:
  /**
   * Constructor.
   */
  PenaltyOperator(Flag operatorType_,
		  ElementLevelSet *elLS_,
		  double factor_,
		  bool penaltyCoeffFlag_,
		  FiniteElemSpace *rowFESpace_,
		  FiniteElemSpace *colFESpace_ = NULL)
    : Operator(operatorType_, rowFESpace_, colFESpace_),
      elLS(elLS_),
      elStatus(ElementLevelSet::LEVEL_SET_UNDEFINED),
      factor(factor_),
      penaltyCoeffFlag(penaltyCoeffFlag_),
      surfaceOp(NULL),
      dim(getRowFESpace()->getMesh()->getDim()),
      degree(getRowFESpace()->getBasisFcts()->getDegree())
  {
    TEST_EXIT(elLS->getLevelSetFct() && elLS->getMesh())
      ("ElementLevelSet not initialized!\n");

    tempCoords = NEW VectorOfFixVecs<DimVec<double> >(dim, dim, NO_INIT);
  };

  /**
   * Destructor.
   */
  ~PenaltyOperator()
  {
    if (surfaceOp) {
      DELETE surfaceOp;
    }
/*     if (dim == 3) { */
      DELETE tempCoords;
/*     } */
  };

  /** 
   * Calculates the element matrix for this ElInfo and adds it multiplied by
   * factor to userMat.
   * In addition to \ref Operator::getElementMatrix(), it distinguishes
   * between elements divided by the boundary (level set zero intersects the 
   * element) and elements lying completely inside or outside the integration 
   * domain.
   */
  void getElementMatrix(const ElInfo *elInfo, 
76
			ElementMatrix& userMat, 
77
78
79
80
81
82
83
84
85
86
			double factor = 1.0);

  /** 
   * Calculates the element vector for this ElInfo and adds it multiplied by
   * factor to userVec.
   * Similar to \ref getElementMatrix it distinguishes between elements 
   * divided by the boundary and elements lying completely inside or outside 
   * the integration domain.
   */
  void getElementVector(const ElInfo *elInfo, 
87
			ElementVector& userVec, 
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
			double factor = 1.0);

protected:

  /**
   * Calculate coefficient of the penalty term.
   */
  double getPenaltyCoeff(const ElInfo *elInfo);

 protected:
  /**
   * Holds level set function and functionalities for intersection point
   * calculation.
   */
  ElementLevelSet *elLS;

  /**
   * Indicator for the position of an element. It shows whether an element 
   * lies inside the integration domain, outside of the integration domain 
   * or on the boundary.
   */
  int elStatus;

  /** 
   * Coefficient of the penalty term is 1/eps. eps depends on the size of
   * element and on factor.
   */
  double eps;

  /**
   * Factor needed for eps.
   */
  double factor;

  /**
   * Flag to indicate whether to use penalty coefficient.
   *   true - use penalty coefficient
   *   false - do not use
   */
  bool penaltyCoeffFlag;

  /**
   * Surface operator for surface integration
   */
  SurfaceOperator *surfaceOp;

  /**
   * Problem Dimension.
   */
  int dim;
 
  /**
   * Degree of basis functions.
   */
  int degree;
 
  /**
   * Variables used for calculation.
   */
  VectorOfFixVecs<DimVec<double> > *tempCoords;
};

#endif  // AMDIS_PENALTYOPERATOR_H