DOFVector.hh 31.3 KB
Newer Older
1
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
2
#include <algorithm>
3
#include <math.h>
4

5 6 7 8 9
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/utility/tag.hpp>
#include <boost/numeric/mtl/utility/category.hpp>
#include <boost/numeric/linear_algebra/identity.hpp>

10 11 12 13 14 15 16 17 18 19 20 21
#include "FixVec.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "Error.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "AbstractFunction.h"
#include "BoundaryManager.h"
#include "Assembler.h"
22
#include "OpenMP.h"
23
#include "Operator.h"
24
#include "Parameters.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
25
#include "Traverse.h"
26

27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
// Defining the interface for MTL4
namespace mtl {

  // Let MTL4 know that DOFVector it is a column vector
  namespace traits {
    template <typename T>
    struct category< AMDiS::DOFVector<T> > 
    {
      typedef tag::dense_col_vector type;
    };
  }

  namespace ashape {
    template <typename T>
    struct ashape< AMDiS::DOFVector<T> > 
    {
      typedef cvec<typename ashape<T>::type> type;
    };
  }

  // Modelling Collection and MutableCollection
  template <typename T>
  struct Collection< AMDiS::DOFVector<T> > 
  {
    typedef T               value_type;
    typedef const T&        const_reference;
    typedef std::size_t     size_type;
  };

  template <typename T>
  struct MutableCollection< AMDiS::DOFVector<T> > 
    : public Collection< AMDiS::DOFVector<T> > 
  {
    typedef T&              reference;
  };


} // namespace mtl



68 69 70
namespace AMDiS {

  template<typename T>
71
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
72 73
    : feSpace(f),
      name(n),
74
      elementVector(f->getBasisFcts()->getNumber()),
Thomas Witkowski's avatar
Thomas Witkowski committed
75
      boundaryManager(NULL)
76
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
77
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
78
    dim = feSpace->getMesh()->getDim();
79

80
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
81
    localVectors.resize(omp_get_overall_max_threads());
82
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    grdTmp.resize(omp_get_overall_max_threads());
84
    D2Phis.resize(omp_get_overall_max_threads());
85

86
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
87 88
      localIndices[i] = new DegreeOfFreedom[this->nBasFcts];
      localVectors[i] = new T[this->nBasFcts];
89 90 91
      grdPhis[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      grdTmp[i] = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
      D2Phis[i] = new DimMat<double>(dim, NO_INIT);
92 93
    }
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
94
  
95
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
96 97 98
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
99 100
      delete [] localIndices[i];
      delete [] localVectors[i];
101 102 103
      delete grdPhis[i];
      delete grdTmp[i];
      delete D2Phis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
104 105 106 107
    }
  }
  
  template<typename T>
108
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
109
    : DOFVectorBase<T>(f, n),
110
      coarsenOperation(NO_OPERATION)
111 112 113 114 115
  {
    init(f, n);
  } 

  template<typename T>
116
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
117 118 119
  {
    name = n;
    feSpace = f;
Thomas Witkowski's avatar
Thomas Witkowski committed
120
    if (feSpace && feSpace->getAdmin())
121
      (feSpace->getAdmin())->addDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
122

123
    this->boundaryManager = new BoundaryManager(f);
124 125 126 127 128
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
129
    if (feSpace && feSpace->getAdmin())
130 131
      (feSpace->getAdmin())->removeDOFIndexed(this);

Thomas Witkowski's avatar
Thomas Witkowski committed
132
    if (this->boundaryManager)
133
      delete this->boundaryManager;
134 135

    vec.clear();
136 137 138 139 140 141 142 143 144 145 146
  }

  template<typename T>
  DOFVector<T> * DOFVector<T>::traverseVector = NULL;

  template<typename T>
  FastQuadrature *DOFVector<T>::quad_fast = NULL;

  template<typename T>
  double DOFVector<T>::norm = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
147 148 149 150
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
151
					  ElInfo *elInfo,
Thomas Witkowski's avatar
Thomas Witkowski committed
152 153 154 155
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

156 157 158 159
    Vector<DegreeOfFreedom> indices(nBasFcts);
    feSpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						feSpace->getAdmin(),
						&indices);
Thomas Witkowski's avatar
Thomas Witkowski committed
160

161
    for (DegreeOfFreedom i = 0; i < nBasFcts; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
162 163 164
      BoundaryCondition *condition = 
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

Thomas Witkowski's avatar
Thomas Witkowski committed
165
      if (!(condition && condition->isDirichlet())) {
166
	DegreeOfFreedom irow = indices[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
167 168 169 170 171 172

	if (add) {
	  (*this)[irow] += factor * elVec[i];
	} else {  
	  (*this)[irow] = factor * elVec[i];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
173 174 175 176
      }
    }
  }

177 178 179 180 181
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
182 183
    checkFeSpace(feSpace, vec);

184
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
185 186
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
188 189 190 191 192 193 194 195 196 197
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

Thomas Witkowski's avatar
Thomas Witkowski committed
198
    checkFeSpace(feSpace, vec);
199 200
    
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
201 202
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
203
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
204 205 206 207 208 209 210 211
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
212
    FUNCNAME("DOFVector<T>::asum()");
213

Thomas Witkowski's avatar
Thomas Witkowski committed
214
    checkFeSpace(feSpace, vec);
215

216
    double nrm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
217 218
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
219
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
220 221 222 223 224 225 226 227
      nrm += abs(*vecIterator);

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
228
    FUNCNAME("DOFVector<T>::sum()");
229

Thomas Witkowski's avatar
Thomas Witkowski committed
230
    checkFeSpace(feSpace, vec);
231

232
    double nrm = 0.0;    
Thomas Witkowski's avatar
Thomas Witkowski committed
233 234
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), 
			 USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
235
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
236 237 238 239 240 241 242 243
      nrm += *vecIterator;

    return(nrm);
  }

  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
244
    FUNCNAME("DOFVector<T>::set()");
245

Thomas Witkowski's avatar
Thomas Witkowski committed
246 247
    checkFeSpace(feSpace, vec);

248
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
249
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
250 251 252 253 254 255 256
      *vecIterator = alpha ; 
  }


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
257
    FUNCNAME("DOFVector<T>::copy()");
Thomas Witkowski's avatar
Thomas Witkowski committed
258 259 260

    checkFeSpace(feSpace, vec);

261
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= feSpace->getAdmin()->getUsedSize())
262
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
263 264
       feSpace->getAdmin()->getUsedSize());
    
265
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
266 267 268 269
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), 
		       USED_DOFS);
    for (vecIterator.reset(), xIterator.reset(); !vecIterator.end();
	 ++vecIterator, ++xIterator)      
270
      *vecIterator = *xIterator; 
271 272 273 274 275 276
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
277 278
    FUNCNAME("DOFVector<T>::min()");
    
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    checkFeSpace(feSpace, vec);
280 281

    T m;    
282
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
283
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
284
      m = std::min(m, *vecIterator);
285

286 287 288 289 290 291
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
292 293
    FUNCNAME("DOFVector<T>::max()");
    
Thomas Witkowski's avatar
Thomas Witkowski committed
294
    checkFeSpace(feSpace, vec);
295

296
    T m;    
297
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
298
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator)
299
      m = std::max(m, *vecIterator);
300

301 302 303
    return m;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
304 305 306 307 308 309
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

310 311 312
  template<typename T>
  void DOFVector<T>::print() const
  {
313
    FUNCNAME("DOFVector<T>::print()");
Thomas Witkowski's avatar
Thomas Witkowski committed
314

315
    const DOFAdmin *admin = NULL;
316
    const char *format;
317

318 319
    if (feSpace) 
      admin = feSpace->getAdmin();
320 321

    MSG("Vec `%s':\n", name.c_str());
322
    int j = 0;
323 324 325 326 327 328 329 330 331
    if (admin) {
      if (admin->getUsedSize() > 100)
	format = "%s(%3d,%10.5le)";
      else if (admin->getUsedSize() > 10)
	format = "%s(%2d,%10.5le)";
      else
	format = "%s(%1d,%10.5le)";

      Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
332
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
333
	if ((j % 3) == 0) {
334 335
	  if (j) 
	    Msg::print("\n");
336
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
337
	} else {
338
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
339
	}
340
	j++;
341
      }
342
      Msg::print("\n");
343
    } else {
344 345
	MSG("no DOFAdmin, print whole vector.\n");
    
346 347 348 349 350 351
	for (int i = 0; i < static_cast<int>( vec.size()); i++) {
	  if ((j % 3) == 0) {
	    if (j) 
	      Msg::print("\n");
	    MSG("(%d,%10.5e)",i,vec[i]);
	  } else {
352
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
353
	  }
354 355 356 357 358 359 360
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

361 362 363 364 365 366 367 368 369
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
370 371 372 373 374 375 376 377 378

  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFESpace()->getBasisFcts());
    int numberOfBasFcts = phi->getNumber();
    T val = 0.0;

379
    for (int i = 0; i < numberOfBasFcts; i++)
380 381 382 383 384 385 386 387
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
388
    FUNCNAME("interpol");
389

390 391 392
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
393

394 395 396 397 398
    if (!this->getFESpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
399

400 401 402 403 404
    if (!(this->getFESpace()->getAdmin())) {
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
	  this->getFESpace()->getName().c_str());
      return;
    }
405
  
406 407 408 409 410
    if (!(this->getFESpace()->getBasisFcts())) {
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
411

412 413 414 415 416
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
417 418 419

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
420 421 422 423 424 425 426
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this->getFESpace()->getMesh(), -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
427 428 429 430 431 432 433
  }

  template<typename T>
  int DOFVector<T>::interpolFct(ElInfo* elinfo)
  {
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
Thomas Witkowski's avatar
Thomas Witkowski committed
434 435 436 437 438 439 440
    const DegreeOfFreedom *dof =  basFct->getLocalIndices(const_cast<Element *>(elinfo->getElement()),
							  admin, NULL);
    const T *inter_val = const_cast<BasisFunction*>(basFct)->interpol(elinfo, 
								      0, 
								      NULL,
								      traverseVector->interFct, 
								      NULL);
441
    int number = basFct->getNumber();
442
    for (int i = 0; i < number; i++)
443 444 445 446 447 448 449 450
      (*traverseVector)[dof[i]] = inter_val[i];

    return 0;
  }

  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    FUNCNAME("DOFVector::Int()");
452 453 454
  
    Mesh* mesh = feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
455 456 457 458
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *q, INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   Int_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::Int_fct(ElInfo *elinfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
476 477 478 479
    double det = elinfo->getDet();
    const T* uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int numPoints = quad_fast->getNumPoints();
    double normT = 0.0;
480

Thomas Witkowski's avatar
Thomas Witkowski committed
481 482
    for (int iq = 0; iq < numPoints; iq++)
      normT += quad_fast->getWeight(iq) * (uh_vec[iq]);
483

Thomas Witkowski's avatar
Thomas Witkowski committed
484
    norm += det * normT;
485 486 487 488 489 490 491

    return 0;
  }

  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
492
    FUNCNAME("DOFVector::L1Norm()");
493
  
Thomas Witkowski's avatar
Thomas Witkowski committed
494 495 496 497
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
498 499 500 501

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),*q,INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
502
    Mesh* mesh = feSpace->getMesh();
503 504 505 506 507 508 509 510 511 512 513 514 515

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   L1Norm_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L1Norm_fct(ElInfo *elinfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
516 517 518 519
    double det = elinfo->getDet();
    const T* uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int numPoints = quad_fast->getNumPoints();
    double normT = 0.0;
520

Thomas Witkowski's avatar
Thomas Witkowski committed
521 522
    for (int iq = 0; iq < numPoints; iq++)
      normT += quad_fast->getWeight(iq) * abs(uh_vec[iq]);
523

Thomas Witkowski's avatar
Thomas Witkowski committed
524
    norm += det * normT;
525 526 527 528 529 530 531 532

    return 0;
  }


  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
533
    FUNCNAME("DOFVector::L2NormSquare()");
534

Thomas Witkowski's avatar
Thomas Witkowski committed
535 536 537 538
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
539 540 541 542 543 544 545

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, 
						      INIT_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
546
    Mesh* mesh = feSpace->getMesh();
Thomas Witkowski's avatar
Thomas Witkowski committed
547 548
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
		   L2Norm_fct);
549 550 551 552 553 554 555

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L2Norm_fct(ElInfo *elinfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
556 557 558 559
    double det = elinfo->getDet();
    const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int numPoints = quad_fast->getNumPoints();
    double normT = 0.0;
560

Thomas Witkowski's avatar
Thomas Witkowski committed
561 562
    for (int iq = 0; iq < numPoints; iq++)
      normT += quad_fast->getWeight(iq) * sqr(uh_vec[iq]);
563

Thomas Witkowski's avatar
Thomas Witkowski committed
564
    norm += det * normT;
565 566 567 568 569 570 571 572

    return 0;
  }

  template<typename T>
  int DOFVector<T>::H1Norm_fct(ElInfo *elinfo)
  {
    double det = elinfo->getDet();
Thomas Witkowski's avatar
Thomas Witkowski committed
573 574
    const WorldVector<T> *grduh_vec = 
      traverseVector->getGrdAtQPs(elinfo, NULL, quad_fast, NULL);
575 576
    int dimOfWorld = Global::getGeo(WORLD);
    int numPoints = quad_fast->getNumPoints();
Thomas Witkowski's avatar
Thomas Witkowski committed
577
    double normT = 0.0;
578

Thomas Witkowski's avatar
Thomas Witkowski committed
579 580 581 582 583 584 585 586
    for (int iq = 0; iq < numPoints; iq++) {
      double norm2 = 0.0;
      for (int j = 0; j < dimOfWorld; j++)
	norm2 += sqr(grduh_vec[iq][j]);
      
      normT += quad_fast->getWeight(iq) * norm2;
    }
    norm += det * normT; 
587 588 589 590 591 592 593 594

    return 0;
  }


  template<typename T>
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
595
    FUNCNAME("DOFVector::H1NormSquare()");
596

Thomas Witkowski's avatar
Thomas Witkowski committed
597 598 599 600 601 602 603
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree() - 2;
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, INIT_GRD_PHI);
604 605

    norm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
606 607
    traverseVector = const_cast<DOFVector<T>*>(this); 
    Mesh *mesh = feSpace->getMesh();
608 609 610 611 612 613 614 615 616 617

    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
		   Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA,
		   H1Norm_fct);

    return norm;
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
618
					std::vector<DegreeOfFreedom> &newDOF)
619
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
620 621
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
622
	vec[newDOF[i]] = vec[i];
623 624 625 626 627 628
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
629 630
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
631
	 op != this->operators.end(); ++op)
632
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
633

634 635 636
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
637 638 639 640 641
  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
642
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
643 644 645
      (*it)->finishAssembling();
  }

646
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
647
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
648
  {
649
    feSpace = rhs.feSpace;
650
    DOFVectorBase<T>::feSpace = rhs.feSpace;
651
    this->nBasFcts = rhs.nBasFcts;
652
    vec = rhs.vec;
653
    this->elementVector.change_dim(this->nBasFcts);
654 655 656 657
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
658
    if (rhs.boundaryManager) {
659 660 661
      if (this->boundaryManager) 
	delete this->boundaryManager; 

662
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
663
    } else {
664
      this->boundaryManager = NULL;
665 666
    }

667 668 669 670 671 672
    return *this;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal)
  {
673
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, T scal)");
674

675 676
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
677

Thomas Witkowski's avatar
Thomas Witkowski committed
678 679
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
680
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
681
      (*vecIterator) *= scal; 
682

683 684 685 686 687 688 689
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
690 691 692
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
693
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
694 695
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
696 697
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
698 699
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
700 701
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
702 703
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
704 705
      *xIterator += *yIterator; 

706 707 708 709 710 711
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y)
  {
712
    FUNCNAME("DOFVector<T>::operator-=(DOFVector<T>& x, const DOFVector<T>& y)");
713

714
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
715
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
716 717
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
718 719
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
720 721
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
722 723
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
724 725
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
726
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
727

728 729 730 731 732 733
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
734 735 736
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
737
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
738 739
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
740 741
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
742 743
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
744 745
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
746 747
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
748
      *xIterator *= *yIterator; 
749 750 751 752 753 754 755 756

    return x;
  }

  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");
757
    const DOFAdmin *admin = NULL;
758

759
    TEST_EXIT(x.getFESpace() && y.getFESpace())
760
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
761
    TEST_EXIT((admin = x.getFESpace()->getAdmin()) && (admin == y.getFESpace()->getAdmin()))
762 763
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
764 765
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
766
    T dot = 0;
767 768 769

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
770
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
771
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
772
      dot += (*xIterator) * (*yIterator);
773

Thomas Witkowski's avatar
Thomas Witkowski committed
774
    return dot;
775 776 777
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
778 779
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
780
  {
781 782 783
    // Unfortunately needed
    using namespace mtl;

784 785
    FUNCNAME("DOFVector<T>::mv()");

Thomas Witkowski's avatar
Thomas Witkowski committed
786 787
    TEST_EXIT(a.getRowFESpace() && a.getColFESpace() && 
	      x.getFESpace() && result.getFESpace())
788 789
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
790 791 792 793 794 795 796 797 798

    const DOFAdmin *rowAdmin = a.getRowFESpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFESpace()->getAdmin();

    TEST_EXIT((rowAdmin && colAdmin) && 
	      (((transpose == NoTranspose) && (colAdmin == x.getFESpace()->getAdmin()) && 
		(rowAdmin == result.getFESpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFESpace()->getAdmin()) && 
		(colAdmin == result.getFESpace()->getAdmin()))))
799 800 801 802 803
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());


804 805 806 807 808
    if (transpose == NoTranspose)
      mult(a.getBaseMatrix(), x, result);
    else if (transpose == Transpose)
      mult(trans(const_cast<DOFMatrix::base_matrix_type&>(a.getBaseMatrix())), x, result);
    else 
809 810 811 812
      ERROR_EXIT("transpose = %d\n", transpose);
  }

  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
813
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
814 815 816
  {
    FUNCNAME("DOFVector<T>::axpy()");

817
    TEST_EXIT(x.getFESpace() && y.getFESpace())
818 819
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

820 821 822
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT((admin) && (admin == y.getFESpace()->getAdmin()))
823 824
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
825
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
826
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
827 828
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
829
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848
       admin->getUsedSize());

      // This is the old implementation of the mv-multiplication. It have been changed
      // because of the OpenMP-parallelization:
//     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
//     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
//     for(xIterator.reset(), yIterator.reset();
// 	!xIterator.end();
// 	++xIterator, ++yIterator) 
//       {  
// 	*yIterator += alpha * (*xIterator); 
//       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
849 850
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
851
	  y[i] += alpha * x[i];
852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882
  }

  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v)
  {
    static DOFVector<T> result;
    return mult(d, v, result);
  }

  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T> &v1 , const DOFVector<T> &v2)
  {
    static DOFVector<T> result;
    return add(v1, v2, result);
  }


  template<typename T>
  void xpay(double alpha, 
	    const DOFVector<T>& x, 
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::xpay()");

883
    TEST_EXIT(x.getFESpace() && y.getFESpace())
884 885
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

886 887 888
    const DOFAdmin *admin = x.getFESpace()->getAdmin();

    TEST_EXIT(admin && (admin == y.getFESpace()->getAdmin()))
889 890
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
891
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
892
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
893 894
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
895
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914
       admin->getUsedSize());

    // This is the old implementation of the mv-multiplication. It have been changed
    // because of the OpenMP-parallelization:
    //     typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
    //     typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
    //     for(xIterator.reset(), yIterator.reset();
    // 	!xIterator.end();
    // 	++xIterator, ++yIterator) 
    //       {  
    // 	*yIterator = alpha *(*yIterator)+ (*xIterator); 
    //       };

      int i;
      int maxI = y.getSize();

#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 25000) default(shared) private(i)
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
915 916
      for (i = 0; i < maxI; i++)
	if (!admin->isDOFFree(i))
917
	  y[i] = alpha * y[i] + x[i];
918 919 920 921 922 923 924 925 926
  }

  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
927 928 929
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
930 931 932 933 934 935 936 937 938 939 940

    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
941 942 943 944
    for (