RobinBC.cc 8.62 KB
Newer Older
1
2
3
4
5
#include "RobinBC.h"
#include "Estimator.h"
#include "Assembler.h"
#include "DOFVector.h"
#include "DOFMatrix.h"
6
#include "OpenMP.h"
7
8
9
10
11
12
13
14
#include "SurfaceOperator.h"
#include <math.h>

namespace AMDiS {

  RobinBC::RobinBC(BoundaryType type,
		   AbstractFunction<double, WorldVector<double> > *j,
		   AbstractFunction<double, WorldVector<double> > *alpha,
15
16
17
		   FiniteElemSpace *rowFeSpace_,
		   FiniteElemSpace *colFeSpace_)
    : BoundaryCondition(type, rowFeSpace_, colFeSpace_), 
18
19
20
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
21
    int dim = rowFeSpace->getMesh()->getDim();
22
23
24

    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
25
    coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
26
27

    // for all element sides
28
    for (int i = 0; i < dim + 1; i++) {
29
30
31
      coords[i] = 
	new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
					     DimVec<double>(dim, DEFAULT_VALUE, 0.0));
32
      // for each vertex of the side
33
34
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
35
36
37
38
	(*coords[i])[k][index] = 1.0;
      }
    }

39
    if (j) {
40
      Operator *jOp = new Operator(rowFeSpace);
41
      jOp->addZeroOrderTerm(new CoordsAtQP_ZOT(j));
Thomas Witkowski's avatar
Thomas Witkowski committed
42
      neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
43
    
44
      for (int i = 0; i < dim + 1; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
45
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
46

Thomas Witkowski's avatar
Thomas Witkowski committed
47
      delete jOp;
48
49
    }

50
    if (alpha) {
51
      Operator *alphaOp = new Operator(rowFeSpace, colFeSpace);
52
      alphaOp->addZeroOrderTerm(new CoordsAtQP_ZOT(alpha));
Thomas Witkowski's avatar
Thomas Witkowski committed
53
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
54

55
      for (int i = 0; i < dim + 1; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
56
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
57
      
Thomas Witkowski's avatar
Thomas Witkowski committed
58
      delete alphaOp;
59
60
61
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
62

63
64
65
  RobinBC::RobinBC(BoundaryType type,
		   DOFVectorBase<double> *j,
		   DOFVectorBase<double> *alpha,
66
67
68
		   FiniteElemSpace *rowFeSpace_,
		   FiniteElemSpace *colFeSpace_)
    : BoundaryCondition(type, rowFeSpace_, colFeSpace_), 
69
70
71
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
72
    int dim = rowFeSpace->getMesh()->getDim();
73
74
75

    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
76
    coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
77
78

    // for all element sides
79
    for (int i = 0; i < dim + 1; i++) {
80
81
82
      coords[i] =
	new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE, 
					     DimVec<double>(dim, DEFAULT_VALUE, 0.0));
83
      // for each vertex of the side
84
85
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
86
87
88
89
	(*coords[i])[k][index] = 1.0;
      }
    }

90
    if (j) {
91
      Operator *jOp = new Operator(rowFeSpace);
92
      jOp->addZeroOrderTerm(new VecAtQP_ZOT(j, NULL));
Thomas Witkowski's avatar
Thomas Witkowski committed
93
      neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
94
    
95
      for (int i = 0; i < dim + 1; i++)
96
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);    
97

Thomas Witkowski's avatar
Thomas Witkowski committed
98
      delete jOp;
99
100
    }

101
    if (alpha) {
102
      Operator *alphaOp = new Operator(rowFeSpace, colFeSpace);
103
      alphaOp->addZeroOrderTerm(new VecAtQP_ZOT(alpha, NULL));
Thomas Witkowski's avatar
Thomas Witkowski committed
104
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
105

106
107
      for (int i = 0; i < dim + 1; i++)
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);     
108

Thomas Witkowski's avatar
Thomas Witkowski committed
109
      delete alphaOp;
110
111
112
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
113

114
115
  RobinBC::RobinBC(BoundaryType type,
		   Operator* jOp, Operator* alphaOp,
116
117
118
		   FiniteElemSpace *rowFeSpace_,
		   FiniteElemSpace *colFeSpace_) 
    : BoundaryCondition(type, rowFeSpace_, colFeSpace_), 
119
120
121
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
122
    int dim = rowFeSpace->getMesh()->getDim();
123
124
125

    // create barycentric coords for each vertex of each side
    const Element *refElement = Global::getReferenceElement(dim);
126
    coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
127
128

    // for all element sides
129
130
131
132
133
    for (int i = 0; i < dim + 1; i++) {
      coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
						       DimVec<double>(dim,
								      DEFAULT_VALUE, 
								      0.0));
134
      // for each vertex of the side
135
136
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
137
138
139
140
	(*coords[i])[k][index] = 1.0;
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
    neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
    robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
143

144
145
    for (int i = 0; i < dim + 1; i++) {
      if (jOp)
Thomas Witkowski's avatar
Thomas Witkowski committed
146
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
147
      if (alphaOp)
Thomas Witkowski's avatar
Thomas Witkowski committed
148
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
149
150
151
    }   
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
152

153
  void RobinBC::fillBoundaryCondition(DOFVectorBase<double>* vector, 
154
155
156
157
				      ElInfo* elInfo,
				      const DegreeOfFreedom* dofIndices,
				      const BoundaryType* localBound,
				      int nBasFcts)
158
159
  {
    FUNCNAME("RobinBC::fillBoundaryCondition()");
160
    TEST_EXIT_DBG(vector->getFeSpace() == rowFeSpace)("invalid row fe space\n");
161
162
163

    int dim = elInfo->getMesh()->getDim();

Thomas Witkowski's avatar
Thomas Witkowski committed
164
165
166
    if (neumannOperators)
      for (int i = 0; i < dim + 1; i++)
	if (elInfo->getBoundary(i) == boundaryType)
167
168
169
	  vector->assemble(1.0, elInfo, localBound, (*neumannOperators)[i]);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
170

171
172
  void RobinBC::fillBoundaryCondition(DOFMatrix* matrix,
				      ElInfo* elInfo,
173
				      const DegreeOfFreedom* dofIndices,
174
175
				      const BoundaryType* localBound,
				      int nBasFcts) 
176
  {
177
    if (robinOperators) {
178
179
180
181
      int dim = elInfo->getMesh()->getDim();

      for (int i = 0; i < dim + 1; i++)
	if (elInfo->getBoundary(i) == boundaryType)
182
183
184
185
	  matrix->assemble(-1.0, elInfo, localBound, (*robinOperators)[i]);
    }
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
186

187
188
189
190
191
  double RobinBC::boundResidual(ElInfo *elInfo,
				DOFMatrix *matrix,
				const DOFVectorBase<double> *dv)
  {
    FUNCNAME("RobinBC::fillBoundaryCondition()");
192
193
    TEST_EXIT(matrix->getRowFeSpace() == rowFeSpace)("invalid row fe space\n");
    TEST_EXIT(matrix->getColFeSpace() == colFeSpace)("invalid col fe space\n");
194
195
196
197
198
199
200
201

    int dim = elInfo->getMesh()->getDim();
    DimVec<double>  lambda(dim, NO_INIT);
    double n_A_grdUh, val = 0.0;
    WorldVector<double> normal;
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    double det = elInfo->getDet();
    bool neumannQuad = false;
202
    const BasisFunction *basFcts = dv->getFeSpace()->getBasisFcts();
203

204
    TEST_EXIT(basFcts == rowFeSpace->getBasisFcts())("invalid basFcts\n");
205

206
    ElementVector uhEl(basFcts->getNumber());
207
208
209
210
211
    dv->getLocalVector(elInfo->getElement(), uhEl);

    TEST_EXIT(neumannOperators || robinOperators)
      ("neither neumann nor robin operators set!\n");

212
    if (!robinOperators) {
213
      neumannQuad = true;
214
215
216
217
218
219
    } else {
      if (neumannOperators) {
	if ((*neumannOperators)[0]->getAssembler(omp_get_thread_num())->
	    getZeroOrderAssembler()->getQuadrature()->getNumPoints() > 
	    (*robinOperators)[0]->getAssembler(omp_get_thread_num())->
	    getZeroOrderAssembler()->getQuadrature()->getNumPoints()) 
220
	  neumannQuad = true;
221
222
223
      }
    }

224
    std::vector<Operator*>::iterator op;
225
    for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
Thomas Witkowski's avatar
Thomas Witkowski committed
226
      (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);        
227

228
    for (int face = 0; face < dim + 1; face++) {
229
230
      elInfo->getNormal(face, normal);

231
      Quadrature *quadrature = 
232
	neumannQuad ? 
233
	(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->
234
	getZeroOrderAssembler()->getQuadrature() :
235
	(*robinOperators)[face]->getAssembler(omp_get_thread_num())->
236
237
	getZeroOrderAssembler()->getQuadrature();

238
      if (elInfo->getBoundary(face) == boundaryType) {
Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
	(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->
	  getZeroOrderAssembler()->initElement(elInfo);
241

242
243
244
	int nPoints = quadrature->getNumPoints();
	mtl::dense_vector<double> uhAtQp(nPoints);
	dv->getVecAtQPs(elInfo, quadrature, NULL, uhAtQp);
Thomas Witkowski's avatar
Thomas Witkowski committed
245
	std::vector<double> f(nPoints, 0.0);
246

247
	if (robinOperators)
248
249
	  (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp, 
						 NULL,  NULL, &f[0], 1.0);
250
	
Thomas Witkowski's avatar
Thomas Witkowski committed
251
252
	std::vector<WorldVector<double> > grdUh(nPoints);
	std::vector<WorldVector<double> > A_grdUh(nPoints);
253

Thomas Witkowski's avatar
Thomas Witkowski committed
254
	for (int iq = 0; iq < nPoints; iq++) {
255
256
257
258
	  A_grdUh[iq].set(0.0);	
	  lambda = quadrature->getLambda(iq);
	  basFcts->evalGrdUh(lambda, Lambda, uhEl, &grdUh[iq]);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
259
	
260
	for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
Thomas Witkowski's avatar
Thomas Witkowski committed
261
	  (*op)->weakEvalSecondOrder(grdUh, A_grdUh);		
262

263
	if (neumannOperators)
Thomas Witkowski's avatar
Thomas Witkowski committed
264
	  (*neumannOperators)[face]->getC(elInfo, nPoints, f);
265

Thomas Witkowski's avatar
Thomas Witkowski committed
266
267
	val = 0.0;
	for (int iq = 0; iq < nPoints; iq++) {
268
269
270
271
272
273
274
275
276
277
	  n_A_grdUh = (normal*A_grdUh[iq]) - f[iq]; 
	  val += quadrature->getWeight(iq)*sqr(n_A_grdUh);
	}
      }
    }

    return det * val;
  }

}