Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

RobinBC.cc 8.92 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
#include "RobinBC.h"
#include "Assembler.h"
#include "DOFVector.h"
#include "DOFMatrix.h"
#include "SurfaceOperator.h"
18 19
#include "est/Estimator.h"

20 21 22 23 24 25 26
#include <math.h>

namespace AMDiS {

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
74

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
125

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

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

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

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

164 165
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
166

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
184

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

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

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

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

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

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

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

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

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

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

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

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

256 257
	int nPoints = quadrature->getNumPoints();
	mtl::dense_vector<double> uhAtQp(nPoints);
258 259
	mtl::dense_vector<WorldVector<double> > grdUhAtQp;
	mtl::dense_vector<WorldMatrix<double> > D2UhAtQp;
260
	dv->getVecAtQPs(elInfo, quadrature, NULL, uhAtQp);
261 262
	ElementVector f(nPoints);
	f = 0.0;
263

264
	if (robinOperators)
265
	  (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp, grdUhAtQp,  D2UhAtQp, f, 1.0);
266
	
Thomas Witkowski's avatar
Thomas Witkowski committed
267 268
	std::vector<WorldVector<double> > grdUh(nPoints);
	std::vector<WorldVector<double> > A_grdUh(nPoints);
269

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

279
	if (neumannOperators)
Thomas Witkowski's avatar
Thomas Witkowski committed
280
	  (*neumannOperators)[face]->getC(elInfo, nPoints, f);
281

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

    return det * val;
  }

}