RobinBC.cc 9.06 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
15
16
17
18
19
20
#include "SurfaceOperator.h"
#include <math.h>

namespace AMDiS {

  RobinBC::RobinBC(BoundaryType type,
		   AbstractFunction<double, WorldVector<double> > *j,
		   AbstractFunction<double, WorldVector<double> > *alpha,
		   FiniteElemSpace *rowFESpace_,
		   FiniteElemSpace *colFESpace_)
    : BoundaryCondition(type, rowFESpace_, colFESpace_), 
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
21
    int dim = rowFESpace->getMesh()->getDim();
22
23
24
25
26
27

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

    // for all element sides
28
29
30
31
32
    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));
33
      // for each vertex of the side
34
35
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
36
37
38
39
	(*coords[i])[k][index] = 1.0;
      }
    }

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

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

51
52
53
    if (alpha) {
      Operator *alphaOp = 
	new Operator(Operator::MATRIX_OPERATOR, rowFESpace, colFESpace);
54
      alphaOp->addZeroOrderTerm(new CoordsAtQP_ZOT(alpha));
Thomas Witkowski's avatar
Thomas Witkowski committed
55
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
56

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

  RobinBC::RobinBC(BoundaryType type,
		   DOFVectorBase<double> *j,
		   DOFVectorBase<double> *alpha,
		   FiniteElemSpace *rowFESpace_,
		   FiniteElemSpace *colFESpace_)
    : BoundaryCondition(type, rowFESpace_, colFESpace_), 
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
73
    int dim = rowFESpace->getMesh()->getDim();
74
75
76
77
78
79

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

    // for all element sides
80
81
82
83
84
85
    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)
						       );
86
      // for each vertex of the side
87
88
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
89
90
91
92
	(*coords[i])[k][index] = 1.0;
      }
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
101
      delete jOp;
102
103
    }

104
105
106
    if (alpha) {
      Operator *alphaOp = 
	new Operator(Operator::MATRIX_OPERATOR, rowFESpace, colFESpace);
107
      alphaOp->addZeroOrderTerm(new VecAtQP_ZOT(alpha, NULL));
Thomas Witkowski's avatar
Thomas Witkowski committed
108
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
109

110
111
      for (int i = 0; i < dim + 1; i++)
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);     
112

Thomas Witkowski's avatar
Thomas Witkowski committed
113
      delete alphaOp;
114
115
116
117
118
119
120
121
122
123
124
    }
  }

  RobinBC::RobinBC(BoundaryType type,
		   Operator* jOp, Operator* alphaOp,
		   FiniteElemSpace *rowFESpace_,
		   FiniteElemSpace *colFESpace_) 
    : BoundaryCondition(type, rowFESpace_, colFESpace_), 
      neumannOperators(NULL), 
      robinOperators(NULL)
  {
125
    int dim = rowFESpace->getMesh()->getDim();
126
127
128
129
130
131

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

    // for all element sides
132
133
134
135
136
    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));
137
      // for each vertex of the side
138
139
      for (int k = 0; k < dim; k++) {
	int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
140
141
142
143
	(*coords[i])[k][index] = 1.0;
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
    neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
    robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
146

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

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

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

166
167
168
    if (neumannOperators) {
      for (int i = 0; i < dim + 1; i++) {
	if (elInfo->getBoundary(i) == boundaryType) {
169
170
171
172
173
174
	  vector->assemble(1.0, elInfo, localBound, (*neumannOperators)[i]);
	}
      }
    }
  }

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

      for (int i = 0; i < dim + 1; i++)
	if (elInfo->getBoundary(i) == boundaryType)
186
187
188
189
190
191
192
193
194
	  matrix->assemble(-1.0, elInfo, localBound, (*robinOperators)[i]);
    }
  }
  
  double RobinBC::boundResidual(ElInfo *elInfo,
				DOFMatrix *matrix,
				const DOFVectorBase<double> *dv)
  {
    FUNCNAME("RobinBC::fillBoundaryCondition()");
195
196
    TEST_EXIT(matrix->getRowFESpace() == rowFESpace)("invalid row fe space\n");
    TEST_EXIT(matrix->getColFESpace() == colFESpace)("invalid col fe space\n");
197
198
199
200
201
202
203
204
205
206
207
208
209
210

    int dim = elInfo->getMesh()->getDim();
    int iq;
    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();
    int numPoints;
    bool neumannQuad = false;
    const BasisFunction *basFcts = dv->getFESpace()->getBasisFcts();

    TEST_EXIT(basFcts == rowFESpace->getBasisFcts())("invalid basFcts\n");

211
    double *uhEl = new double[basFcts->getNumber()];
212
213
214
215
216
217

    dv->getLocalVector(elInfo->getElement(), uhEl);

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

218
    if (!robinOperators) {
219
      neumannQuad = true;
220
221
222
223
224
225
    } else {
      if (neumannOperators) {
	if ((*neumannOperators)[0]->getAssembler(omp_get_thread_num())->
	    getZeroOrderAssembler()->getQuadrature()->getNumPoints() > 
	    (*robinOperators)[0]->getAssembler(omp_get_thread_num())->
	    getZeroOrderAssembler()->getQuadrature()->getNumPoints()) 
226
227
228
229
230
231
232
233
234
	  {
	    neumannQuad = true;
	  }
      }
    }

    Quadrature *quadrature = NULL;
    int face;

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

240
    for (face = 0; face < dim + 1; face++) {
241
242
243
244
245

      elInfo->getNormal(face, normal);

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

      numPoints = quadrature->getNumPoints();

253
      if (elInfo->getBoundary(face) == boundaryType) {
254
	(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->getZeroOrderAssembler()->
255
256
257
258
259
260
261
	  initElement(elInfo);

	const double *uhAtQp = dv->getVecAtQPs(elInfo,
					       quadrature,
					       NULL,
					       NULL);

262
	std::vector<double> f(numPoints, 0.0);
263

264
	if (robinOperators)
265
266
267
268
	  (*robinOperators)[face]->evalZeroOrder(numPoints, 
						 uhAtQp,
						 NULL,
						 NULL,
269
						 &f[0],
270
						 1.0);
271
	
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
	WorldVector<double> *grdUh = new WorldVector<double>[numPoints];
	WorldVector<double> *A_grdUh = new WorldVector<double>[numPoints];
274
275
276
277
278
279
280

	for (iq = 0; iq < numPoints; iq++) {
	  A_grdUh[iq].set(0.0);	
	  lambda = quadrature->getLambda(iq);
	  basFcts->evalGrdUh(lambda, Lambda, uhEl, &grdUh[iq]);
	}

281
	for (op = matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op) {
282
283
284
285
286
287
	  (*op)->weakEvalSecondOrder(numPoints, 
				     grdUh, 
				     A_grdUh);
	}


288
	if (neumannOperators)
289
290
291
292
293
294
295
	  (*neumannOperators)[face]->getC(elInfo, numPoints, f);

	for (val = iq = 0; iq < numPoints; iq++) {
	  n_A_grdUh = (normal*A_grdUh[iq]) - f[iq]; 
	  val += quadrature->getWeight(iq)*sqr(n_A_grdUh);
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
	delete [] grdUh;
	delete [] A_grdUh;
298
299
300
      }
    }

301
    delete [] uhEl;
302
303
304
305
306

    return det * val;
  }

}