RobinBC.cc 9.28 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
21
22
23
24
25
26
27
28
29
30
31
#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)
		  //matrix(dofMatrix)
  {
    int i, k, dim = rowFESpace->getMesh()->getDim();

    // 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
    for(i=0; i < dim+1; i++) {
      coords[i] = 
Thomas Witkowski's avatar
Thomas Witkowski committed
32
	new VectorOfFixVecs<DimVec<double> >(dim,
33
34
35
36
37
38
39
40
41
42
43
44
45
46
					     dim,
					     DEFAULT_VALUE,
					     DimVec<double>(dim, 
							    DEFAULT_VALUE, 
							    0.0));
      // for each vertex of the side
      for(k=0; k < dim; k++) {
	int index = 
	  refElement->getVertexOfPosition(INDEX_OF_DIM(dim-1, dim), i, k);
	(*coords[i])[k][index] = 1.0;
      }
    }

    if(j) {
Thomas Witkowski's avatar
Thomas Witkowski committed
47
      Operator *jOp = new Operator(Operator::VECTOR_OPERATOR, rowFESpace);
48
      jOp->addZeroOrderTerm(new CoordsAtQP_ZOT(j));
Thomas Witkowski's avatar
Thomas Witkowski committed
49
      neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
50
51
    
      for(i=0; i < dim+1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
52
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
53
54
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
55
      delete jOp;
56
57
58
    }

    if(alpha) {
Thomas Witkowski's avatar
Thomas Witkowski committed
59
      Operator *alphaOp = new Operator(Operator::MATRIX_OPERATOR, 
60
61
				       rowFESpace, colFESpace);
      alphaOp->addZeroOrderTerm(new CoordsAtQP_ZOT(alpha));
Thomas Witkowski's avatar
Thomas Witkowski committed
62
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
63
64

      for(i=0; i < dim + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
65
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
66
67
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
68
      delete alphaOp;
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
    }
  }

  RobinBC::RobinBC(BoundaryType type,
		   DOFVectorBase<double> *j,
		   DOFVectorBase<double> *alpha,
		   FiniteElemSpace *rowFESpace_,
		   FiniteElemSpace *colFESpace_)
    : BoundaryCondition(type, rowFESpace_, colFESpace_), 
      neumannOperators(NULL), 
      robinOperators(NULL)
		  //matrix(dofMatrix)
  {
    int i, k, dim = rowFESpace->getMesh()->getDim();

    // 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
    for(i=0; i < dim+1; i++) {
      coords[i] = 
Thomas Witkowski's avatar
Thomas Witkowski committed
92
	new VectorOfFixVecs<DimVec<double> >(dim, 
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
					     dim, 
					     DEFAULT_VALUE, 
					     DimVec<double>(dim, 
							    DEFAULT_VALUE, 
							    0.0)
					     );
      // for each vertex of the side
      for(k=0; k < dim; k++) {
	int index = 
	  refElement->getVertexOfPosition(INDEX_OF_DIM(dim-1, dim), i, k);
	(*coords[i])[k][index] = 1.0;
      }
    }

    if(j) {
Thomas Witkowski's avatar
Thomas Witkowski committed
108
      Operator *jOp = new Operator(Operator::VECTOR_OPERATOR, rowFESpace);
109
      jOp->addZeroOrderTerm(new VecAtQP_ZOT(j, NULL));
Thomas Witkowski's avatar
Thomas Witkowski committed
110
      neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
111
112
    
      for(i=0; i < dim+1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
113
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
114
115
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
116
      delete jOp;
117
118
119
    }

    if(alpha) {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      Operator *alphaOp = new Operator(Operator::MATRIX_OPERATOR, 
121
122
				       rowFESpace, colFESpace);
      alphaOp->addZeroOrderTerm(new VecAtQP_ZOT(alpha, NULL));
Thomas Witkowski's avatar
Thomas Witkowski committed
123
      robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
124
125

      for(i=0; i < dim + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
126
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
127
128
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
129
      delete alphaOp;
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
    }
  }

  RobinBC::RobinBC(BoundaryType type,
		   Operator* jOp, Operator* alphaOp,
		   FiniteElemSpace *rowFESpace_,
		   FiniteElemSpace *colFESpace_) 
    : BoundaryCondition(type, rowFESpace_, colFESpace_), 
      neumannOperators(NULL), 
      robinOperators(NULL)
		  //matrix(dofMatrix)
  {
    int i, k, dim = rowFESpace->getMesh()->getDim();

    // 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
    for(i=0; i < dim+1; i++) {
      coords[i] = 
Thomas Witkowski's avatar
Thomas Witkowski committed
152
	new VectorOfFixVecs<DimVec<double> >(dim,
153
154
155
156
157
158
159
160
161
162
163
164
165
					     dim, 
					     DEFAULT_VALUE,
					     DimVec<double>(dim,
							    DEFAULT_VALUE, 
							    0.0));
      // for each vertex of the side
      for(k=0; k < dim; k++) {
	int index = 
	  refElement->getVertexOfPosition(INDEX_OF_DIM(dim-1, dim), i, k);
	(*coords[i])[k][index] = 1.0;
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
166
167
    neumannOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
    robinOperators = new DimVec<SurfaceOperator*>(dim, NO_INIT);
168
169
170

    for(i=0; i < dim + 1; i++) {
      if(jOp)
Thomas Witkowski's avatar
Thomas Witkowski committed
171
	(*neumannOperators)[i] = new SurfaceOperator(jOp, *coords[i]);
172
      if(alphaOp)
Thomas Witkowski's avatar
Thomas Witkowski committed
173
	(*robinOperators)[i] = new SurfaceOperator(alphaOp, *coords[i]);
174
175
176
177
    }   
  }

  void RobinBC::fillBoundaryCondition(DOFVectorBase<double>* vector, 
178
179
180
181
				      ElInfo* elInfo,
				      const DegreeOfFreedom* dofIndices,
				      const BoundaryType* localBound,
				      int nBasFcts)
182
183
  {
    FUNCNAME("RobinBC::fillBoundaryCondition()");
184
    TEST_EXIT_DBG(vector->getFESpace() == rowFESpace)("invalid row fe space\n");
185
186
187

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

188
189
190
    if (neumannOperators) {
      for (int i = 0; i < dim + 1; i++) {
	if (elInfo->getBoundary(i) == boundaryType) {
191
192
193
194
195
196
	  vector->assemble(1.0, elInfo, localBound, (*neumannOperators)[i]);
	}
      }
    }
  }

197
198
  void RobinBC::fillBoundaryCondition(DOFMatrix* matrix,
				      ElInfo* elInfo,
199
				      const DegreeOfFreedom* dofIndices,
200
201
				      const BoundaryType* localBound,
				      int nBasFcts) 
202
203
204
  {
    int dim = elInfo->getMesh()->getDim();

205
206
207
    if (robinOperators) {
      for (int i = 0; i < dim + 1; i++) {
	if (elInfo->getBoundary(i) == boundaryType) {
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
	  matrix->assemble(-1.0, elInfo, localBound, (*robinOperators)[i]);
	}
      }
    }
  }
  
  double RobinBC::boundResidual(ElInfo *elInfo,
				DOFMatrix *matrix,
				const DOFVectorBase<double> *dv)
  {
    FUNCNAME("RobinBC::fillBoundaryCondition()");
    TEST_EXIT(matrix->getRowFESpace() == rowFESpace)
      ("invalid row fe space\n");
    TEST_EXIT(matrix->getColFESpace() == colFESpace)
      ("invalid col fe space\n");

    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");

    double *uhEl = GET_MEMORY(double, basFcts->getNumber());

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

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

    if(!robinOperators) 
      neumannQuad = true;
    else {
      if(neumannOperators) {
252
	if((*neumannOperators)[0]->getAssembler(omp_get_thread_num())->
253
	   getZeroOrderAssembler()->getQuadrature()->getNumPoints() > 
254
	   (*robinOperators)[0]->getAssembler(omp_get_thread_num())->
255
256
257
258
259
260
261
262
263
264
265
	   getZeroOrderAssembler()->getQuadrature()->getNumPoints()) 
	  {
	    neumannQuad = true;
	  }
      }
    }

    Quadrature *quadrature = NULL;

    int face;

266
    std::vector<Operator*>::iterator op;
267
    for(op=matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op) {
268
      (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
269
270
271
272
273
274
275
276
    }

    for(face = 0; face < dim+1; face++) {

      elInfo->getNormal(face, normal);

      quadrature = 
	neumannQuad ? 
277
	(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->
278
	getZeroOrderAssembler()->getQuadrature() :
279
	(*robinOperators)[face]->getAssembler(omp_get_thread_num())->
280
281
282
283
284
285
	getZeroOrderAssembler()->getQuadrature();

      numPoints = quadrature->getNumPoints();

      if(elInfo->getBoundary(face) == boundaryType) {

286
	(*neumannOperators)[face]->getAssembler(omp_get_thread_num())->getZeroOrderAssembler()->
287
288
289
290
291
292
293
294
	  initElement(elInfo);

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

295
	std::vector<double> f(numPoints, 0.0);
296
297
298
299
300
301

	if(robinOperators) {
	  (*robinOperators)[face]->evalZeroOrder(numPoints, 
						 uhAtQp,
						 NULL,
						 NULL,
302
						 &f[0],
303
304
305
						 1.0);
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
	WorldVector<double> *grdUh = new WorldVector<double>[numPoints];
	WorldVector<double> *A_grdUh = new WorldVector<double>[numPoints];
308
309
310
311
312
313
314
315
316
317
318
319
320
321

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

	for(op=matrix->getOperatorsBegin(); op != matrix->getOperatorsEnd(); ++op) {
	  (*op)->weakEvalSecondOrder(numPoints, 
				     grdUh, 
				     A_grdUh);
	}


322
	if (neumannOperators)
323
324
325
326
327
328
329
	  (*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
330
331
	delete [] grdUh;
	delete [] A_grdUh;
332
333
334
335
336
337
338
339
340
      }
    }

    FREE_MEMORY(uhEl, double, basFcts->getNumber());

    return det * val;
  }

}