RobinBC.cc 8.78 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
17
18
19
20
21
22
23
24
25
#include "RobinBC.h"
#include "Estimator.h"
#include "Assembler.h"
#include "DOFVector.h"
#include "DOFMatrix.h"
#include "SurfaceOperator.h"
#include <math.h>

namespace AMDiS {

  RobinBC::RobinBC(BoundaryType type,
		   AbstractFunction<double, WorldVector<double> > *j,
		   AbstractFunction<double, WorldVector<double> > *alpha,
26
27
28
		   FiniteElemSpace *rowFeSpace_,
		   FiniteElemSpace *colFeSpace_)
    : BoundaryCondition(type, rowFeSpace_, colFeSpace_), 
29
30
31
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
32
    int dim = rowFeSpace->getMesh()->getDim();
33
34
35

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

    // for all element sides
39
    for (int i = 0; i < dim + 1; i++) {
40
41
42
      coords[i] = 
	new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
					     DimVec<double>(dim, DEFAULT_VALUE, 0.0));
43
      // for each vertex of the side
44
45
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
46
47
48
49
	(*coords[i])[k][index] = 1.0;
      }
    }

50
    if (j) {
51
      Operator *jOp = new Operator(rowFeSpace);
52
      jOp->addZeroOrderTerm(new CoordsAtQP_ZOT(j));
Thomas Witkowski's avatar
Thomas Witkowski committed
53
      neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
54
    
55
      for (int i = 0; i < dim + 1; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
56
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
57

Thomas Witkowski's avatar
Thomas Witkowski committed
58
      delete jOp;
59
60
    }

61
    if (alpha) {
62
      Operator *alphaOp = new Operator(rowFeSpace, colFeSpace);
63
      alphaOp->addZeroOrderTerm(new CoordsAtQP_ZOT(alpha));
Thomas Witkowski's avatar
Thomas Witkowski committed
64
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
65

66
      for (int i = 0; i < dim + 1; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
67
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
68
      
Thomas Witkowski's avatar
Thomas Witkowski committed
69
      delete alphaOp;
70
71
72
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
73

74
75
76
  RobinBC::RobinBC(BoundaryType type,
		   DOFVectorBase<double> *j,
		   DOFVectorBase<double> *alpha,
77
78
79
		   FiniteElemSpace *rowFeSpace_,
		   FiniteElemSpace *colFeSpace_)
    : BoundaryCondition(type, rowFeSpace_, colFeSpace_), 
80
81
82
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
83
    int dim = rowFeSpace->getMesh()->getDim();
84
85
86

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

    // for all element sides
90
    for (int i = 0; i < dim + 1; i++) {
91
92
93
      coords[i] =
	new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE, 
					     DimVec<double>(dim, DEFAULT_VALUE, 0.0));
94
      // for each vertex of the side
95
96
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
97
98
99
100
	(*coords[i])[k][index] = 1.0;
      }
    }

101
    if (j) {
102
      Operator *jOp = new Operator(rowFeSpace);
103
      jOp->addZeroOrderTerm(new VecAtQP_ZOT(j, NULL));
Thomas Witkowski's avatar
Thomas Witkowski committed
104
      neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
105
    
106
      for (int i = 0; i < dim + 1; i++)
107
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);    
108

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

112
    if (alpha) {
113
      Operator *alphaOp = new Operator(rowFeSpace, colFeSpace);
114
      alphaOp->addZeroOrderTerm(new VecAtQP_ZOT(alpha, NULL));
Thomas Witkowski's avatar
Thomas Witkowski committed
115
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
116

117
118
      for (int i = 0; i < dim + 1; i++)
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);     
119

Thomas Witkowski's avatar
Thomas Witkowski committed
120
      delete alphaOp;
121
122
123
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
124

125
126
  RobinBC::RobinBC(BoundaryType type,
		   Operator* jOp, Operator* alphaOp,
127
128
129
		   FiniteElemSpace *rowFeSpace_,
		   FiniteElemSpace *colFeSpace_) 
    : BoundaryCondition(type, rowFeSpace_, colFeSpace_), 
130
131
132
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
133
    int dim = rowFeSpace->getMesh()->getDim();
134
135
136

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

    // for all element sides
140
141
142
143
144
    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));
145
      // for each vertex of the side
146
147
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
148
149
150
151
	(*coords[i])[k][index] = 1.0;
      }
    }

152
153
154
    if (jOp) {
      neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
      for (int i = 0; i < dim + 1; i++) 
Thomas Witkowski's avatar
Thomas Witkowski committed
155
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
156
157
158
159
    }
    if (alphaOp) {
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
      for (int i = 0; i < dim + 1; i++) 
Thomas Witkowski's avatar
Thomas Witkowski committed
160
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
161
162
    }

163
164
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
165

166
  void RobinBC::fillBoundaryCondition(DOFVectorBase<double>* vector, 
167
168
169
170
				      ElInfo* elInfo,
				      const DegreeOfFreedom* dofIndices,
				      const BoundaryType* localBound,
				      int nBasFcts)
171
172
  {
    FUNCNAME("RobinBC::fillBoundaryCondition()");
173
    TEST_EXIT_DBG(vector->getFeSpace() == rowFeSpace)("invalid row fe space\n");
174
175
176

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

Thomas Witkowski's avatar
Thomas Witkowski committed
177
178
179
    if (neumannOperators)
      for (int i = 0; i < dim + 1; i++)
	if (elInfo->getBoundary(i) == boundaryType)
180
181
182
	  vector->assemble(1.0, elInfo, localBound, (*neumannOperators)[i]);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
183

184
185
  void RobinBC::fillBoundaryCondition(DOFMatrix* matrix,
				      ElInfo* elInfo,
186
				      const DegreeOfFreedom* dofIndices,
187
188
				      const BoundaryType* localBound,
				      int nBasFcts) 
189
  {
190
    if (robinOperators) {
191
192
193
194
      int dim = elInfo->getMesh()->getDim();

      for (int i = 0; i < dim + 1; i++)
	if (elInfo->getBoundary(i) == boundaryType)
195
196
197
198
	  matrix->assemble(-1.0, elInfo, localBound, (*robinOperators)[i]);
    }
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
199

200
201
202
203
204
  double RobinBC::boundResidual(ElInfo *elInfo,
				DOFMatrix *matrix,
				const DOFVectorBase<double> *dv)
  {
    FUNCNAME("RobinBC::fillBoundaryCondition()");
205
206
    TEST_EXIT(matrix->getRowFeSpace() == rowFeSpace)("invalid row fe space\n");
    TEST_EXIT(matrix->getColFeSpace() == colFeSpace)("invalid col fe space\n");
207
208
209
210
211
212
213
214

    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;
215
    const BasisFunction *basFcts = dv->getFeSpace()->getBasisFcts();
216

217
    TEST_EXIT(basFcts == rowFeSpace->getBasisFcts())("invalid basFcts\n");
218

219
    ElementVector uhEl(basFcts->getNumber());
220
221
222
223
224
    dv->getLocalVector(elInfo->getElement(), uhEl);

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

225
    if (!robinOperators) {
226
      neumannQuad = true;
227
228
    } else {
      if (neumannOperators) {
229
	if ((*neumannOperators)[0]->getAssembler()->
230
	    getZeroOrderAssembler()->getQuadrature()->getNumPoints() > 
231
	    (*robinOperators)[0]->getAssembler()->
232
	    getZeroOrderAssembler()->getQuadrature()->getNumPoints()) 
233
	  neumannQuad = true;
234
235
236
      }
    }

237
    std::vector<Operator*>::iterator op;
238
    for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
239
      (*op)->getAssembler()->initElement(elInfo);        
240

241
    for (int face = 0; face < dim + 1; face++) {
242
243
      elInfo->getNormal(face, normal);

244
      Quadrature *quadrature = 
245
	neumannQuad ? 
246
	(*neumannOperators)[face]->getAssembler()->
247
	getZeroOrderAssembler()->getQuadrature() :
248
	(*robinOperators)[face]->getAssembler()->
249
250
	getZeroOrderAssembler()->getQuadrature();

251
      if (elInfo->getBoundary(face) == boundaryType) {
252
	(*neumannOperators)[face]->getAssembler()->
Thomas Witkowski's avatar
Thomas Witkowski committed
253
	  getZeroOrderAssembler()->initElement(elInfo);
254

255
256
257
	int nPoints = quadrature->getNumPoints();
	mtl::dense_vector<double> uhAtQp(nPoints);
	dv->getVecAtQPs(elInfo, quadrature, NULL, uhAtQp);
258
259
	ElementVector f(nPoints);
	f = 0.0;
260

261
	if (robinOperators)
262
	  (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp, NULL,  NULL, f, 1.0);
263
	
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
	std::vector<WorldVector<double> > grdUh(nPoints);
	std::vector<WorldVector<double> > A_grdUh(nPoints);
266

Thomas Witkowski's avatar
Thomas Witkowski committed
267
	for (int iq = 0; iq < nPoints; iq++) {
268
269
270
271
	  A_grdUh[iq].set(0.0);	
	  lambda = quadrature->getLambda(iq);
	  basFcts->evalGrdUh(lambda, Lambda, uhEl, &grdUh[iq]);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
272
	
273
	for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op)
Thomas Witkowski's avatar
Thomas Witkowski committed
274
	  (*op)->weakEvalSecondOrder(grdUh, A_grdUh);		
275

276
	if (neumannOperators)
Thomas Witkowski's avatar
Thomas Witkowski committed
277
	  (*neumannOperators)[face]->getC(elInfo, nPoints, f);
278

Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
	val = 0.0;
	for (int iq = 0; iq < nPoints; iq++) {
281
282
	  n_A_grdUh = (normal * A_grdUh[iq]) - f[iq]; 
	  val += quadrature->getWeight(iq) * sqr(n_A_grdUh);
283
284
285
286
287
288
289
290
	}
      }
    }

    return det * val;
  }

}