* First part of merge with branches/parallel_assemble

AMDiS/src/Estimator.cc

@@ -100,6 +100,7 @@ namespace AMDiS {
   
     TraverseStack stack;
     ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
+
     while(elInfo) {
       estimateElement(elInfo);
       elInfo = stack.traverseNext(elInfo);

AMDiS/src/Operator.hh

+namespace AMDiS {
+
+  template <typename T>
+  void Operator::addSecondOrderTerm(T *term)
+  {
+    secondOrder[0].push_back(term);
+    term->operat = this;
+
+    for (int i = 1; i < omp_get_max_threads(); i++) {
+      T *newTerm = NEW T(static_cast<const T>(*term));
+      secondOrder[i].push_back(newTerm);
+    }
+  }
+
+  template <typename T>
+  void Operator::addFirstOrderTerm(T *term,
+				   FirstOrderType type)
+  {
+    if (type == GRD_PSI) {
+      firstOrderGrdPsi[0].push_back(term);
+    } else {
+      firstOrderGrdPhi[0].push_back(term);
+    }
+    term->operat = this;
+
+    for (int i = 1; i < omp_get_max_threads(); i++) {
+      T *newTerm = NEW T(static_cast<const T>(*term));
+      if (type == GRD_PSI) {
+	firstOrderGrdPsi[i].push_back(newTerm);
+      } else {
+	firstOrderGrdPhi[i].push_back(newTerm);
+      } 
+    }
+  }
+
+  template <typename T>
+  void Operator::addZeroOrderTerm(T *term)
+  {
+    zeroOrder[0].push_back(term);
+    term->operat = this;
+
+    for (int i = 1; i < omp_get_max_threads(); i++) {
+      T *newTerm = NEW T(static_cast<const T>(*term));
+      zeroOrder[i].push_back(newTerm);
+    }
+  }
+
+
+}

AMDiS/src/PardisoSolver.cc

@@ -125,8 +125,8 @@ namespace AMDiS {
     iparm[7] = 2; // Max numbers of iterative refinement steps
     iparm[9] = 13; // Perturb the pivot elements with 1e-13
     iparm[10] = 1; // Use nonsymmetric permutation and scaling MPS
-    iparm[17] = -1; // Output: Number of nonzeros in the factor LU
-    iparm[18] = -1; // Output: Mflops for LU factorization
+    iparm[17] = 0; // Output: Number of nonzeros in the factor LU
+    iparm[18] = 0; // Output: Mflops for LU factorization
 
     // Maximum number of numerical factorizations
     MKL_INT maxfct = 1; 
@@ -135,7 +135,7 @@ namespace AMDiS {
     MKL_INT mnum = 1;
 
     // Print statistical information in file
-    MKL_INT msglvl = 1;
+    MKL_INT msglvl = 0;
 
     // Error flag
     MKL_INT error = 0;

AMDiS/src/ResidualEstimator.cc

+#include "ResidualEstimator.h"
+#include "Operator.h"
+#include "DOFMatrix.h"
+#include "DOFVector.h"
+#include "Assembler.h"
+#include "Traverse.h"
+#include "Parameters.h"
+
+namespace AMDiS {
+
+  ResidualEstimator::ResidualEstimator(::std::string name, int r) 
+    : Estimator(name, r),
+      C0(1.0), 
+      C1(1.0), 
+      C2(1.0), 
+      C3(1.0)
+  {
+    GET_PARAMETER(0, name + "->C0", "%f", &C0);
+    GET_PARAMETER(0, name + "->C1", "%f", &C1);
+    GET_PARAMETER(0, name + "->C2", "%f", &C2);
+    GET_PARAMETER(0, name + "->C3", "%f", &C3);
+
+    C0 = C0 > 1.e-25 ? sqr(C0) : 0.0;
+    C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
+    C2 = C2 > 1.e-25 ? sqr(C2) : 0.0;
+    C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;
+  }
+
+  void ResidualEstimator::init(double ts)
+  {
+    FUNCNAME("ResidualEstimator::init()");
+
+    timestep = ts;
+
+    mesh = uh[row == -1 ? 0 : row]->getFESpace()->getMesh();
+
+    numSystems = static_cast<int>(uh.size());
+    TEST_EXIT_DBG(numSystems > 0)("no system set\n");
+
+    dim = mesh->getDim();
+    basFcts = GET_MEMORY(const BasisFunction*, numSystems);
+    quadFast = GET_MEMORY(FastQuadrature*, numSystems);
+
+    degree = 0;
+    for (int system = 0; system < numSystems; system++) {
+      basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
+      degree = ::std::max(degree, basFcts[system]->getDegree());
+    }
+
+    degree *= 2;
+
+    quad = Quadrature::provideQuadrature(dim, degree);
+    numPoints = quad->getNumPoints();
+
+    Flag flag = INIT_PHI | INIT_GRD_PHI;
+    if (degree > 2) {
+      flag |= INIT_D2_PHI;
+    }
+
+    for (int system = 0; system < numSystems; system++) {
+      quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system], 
+							       *quad, 
+							       flag);
+    }
+  
+    uhEl = GET_MEMORY(double*, numSystems);
+    uhOldEl = timestep ? GET_MEMORY(double*, numSystems) : NULL;
+
+    for (int system = 0; system < numSystems; system++) {
+      uhEl[system] = GET_MEMORY(double, basFcts[system]->getNumber()); 
+      if (timestep)
+	uhOldEl[system] = GET_MEMORY(double, basFcts[system]->getNumber());
+    }
+
+    uhQP = timestep ? GET_MEMORY(double, numPoints) : NULL;
+    uhOldQP = timestep ? GET_MEMORY(double, numPoints) : NULL;
+
+    riq = GET_MEMORY(double, numPoints);
+
+    TraverseStack stack;
+    ElInfo *elInfo = NULL;
+
+    // clear error indicators and mark elements for jumpRes
+    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
+    while (elInfo) {
+      elInfo->getElement()->setEstimation(0.0, row);
+      elInfo->getElement()->setMark(1);
+      elInfo = stack.traverseNext(elInfo);
+    }
+
+    est_sum = 0.0;
+    est_max = 0.0;
+    est_t_sum = 0.0;
+    est_t_max = 0.0;
+
+    traverseFlag = 
+      Mesh::FILL_NEIGH      |
+      Mesh::FILL_COORDS     |
+      Mesh::FILL_OPP_COORDS |
+      Mesh::FILL_BOUND      |
+      Mesh::FILL_GRD_LAMBDA |
+      Mesh::FILL_DET        |
+      Mesh::CALL_LEAF_EL;
+  }
+
+  void ResidualEstimator::exit(bool output)
+  {
+    FUNCNAME("ResidualEstimator::exit()");
+
+    est_sum = sqrt(est_sum);
+    est_t_sum = sqrt(est_t_sum);
+
+    for (int system = 0; system < numSystems; system++) {
+      FREE_MEMORY(uhEl[system], double, basFcts[system]->getNumber());
+      if (timestep)
+	FREE_MEMORY(uhOldEl[system], double, basFcts[system]->getNumber());    
+    }
+
+    FREE_MEMORY(uhEl, double*, numSystems);
+    if (timestep)
+      FREE_MEMORY(uhOldEl, double*, numSystems);
+
+    if (timestep) {
+      FREE_MEMORY(uhQP, double, numPoints);
+      FREE_MEMORY(uhOldQP, double, numPoints);
+    }
+
+    if (output) {
+      MSG("estimate   = %.8e\n", est_sum);
+      if (C3) {
+	MSG("time estimate   = %.8e\n", est_t_sum);
+      }
+    }
+
+    FREE_MEMORY(riq, double, numPoints);
+    FREE_MEMORY(basFcts, const BasisFunction*, numSystems);
+    FREE_MEMORY(quadFast, FastQuadrature*, numSystems);
+  }
+
+  void ResidualEstimator::estimateElement(ElInfo *elInfo)
+  {
+    FUNCNAME("ResidualEstimator::estimateElement()");
+
+    double est_el = 0.0;
+    double val = 0.0;
+    Element *el, *neigh;
+    int iq;
+
+    TEST_EXIT_DBG(numSystems > 0)("no system set\n");
+
+    ::std::vector<Operator*>::iterator it;
+
+    WorldVector<double> *grdUh_qp = NULL;
+    WorldMatrix<double> *D2uhqp = NULL; 
+
+    el = elInfo->getElement();
+
+    double det = elInfo->getDet();
+    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
+
+    est_el = el->getEstimation(row);
+
+    double h2 = h2_from_det(det, dim);
+
+    for (iq = 0; iq < numPoints; iq++) {
+      riq[iq] = 0.0;
+    }
+
+    for (int system = 0; system < numSystems; system++) {
+
+      if (matrix[system] == NULL) 
+	continue;
+
+      // init assemblers
+      ::std::vector<Operator*>::iterator it;
+
+      for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin();
+	   it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
+	   ++it) {
+	(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, quad);
+      }
+
+      for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
+	   it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd(); 
+	   ++it) {
+	(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, quad);
+      }
+
+      if (timestep) {
+	TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
+	uhOld[system]->getLocalVector(el, uhOldEl[system]);
+  
+	// ===== time and element residuals       
+	if (C0 || C3) {   
+	  uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
+	  uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
+	  
+	  if (C3 && uhOldQP && system == ::std::max(row, 0)) {
+	    for (val = iq = 0; iq < numPoints; iq++) {
+	      double tiq = (uhQP[iq] - uhOldQP[iq]);
+	      val += quad->getWeight(iq) * tiq * tiq;
+	    }
+	    double v = C3 * det * val;
+	    est_t_sum += v;
+	    est_t_max = max(est_t_max, v);
+	  }
+	}
+      }
+      
+      
+      if (C0) {  
+	for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(); 
+	     it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
+	     ++it) {
+	  if ((*it)->zeroOrderTerms() && !uhQP) {
+	    uhQP = GET_MEMORY(double, numPoints);
+	    uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
+	  }
+	  if ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi() 
+	      && !grdUh_qp) {
+	    grdUh_qp = new WorldVector<double>[numPoints];
+	    uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
+	  }
+	  if ((*it)->secondOrderTerms() && !D2uhqp) { 
+	    if (degree > 2) {
+	      D2uhqp = new WorldMatrix<double>[numPoints];
+	      uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
+	    }
+	  }
+	}
+	
+	r(elInfo,
+	  numPoints, 
+	  uhQP,
+	  grdUh_qp,
+	  D2uhqp,
+	  uhOldQP,
+	  NULL,  // grdUhOldQP 
+	  NULL,  // D2UhOldQP
+	  matrix[system], 
+	  fh[system],
+	  quad,
+	  riq);			
+      }     
+    }
+
+    // add integral over r square
+    for (val = iq = 0; iq < numPoints; iq++) {
+      val += quad->getWeight(iq) * riq[iq] * riq[iq];
+    }
+    
+    if (timestep != 0.0 || norm == NO_NORM || norm == L2_NORM)
+      val = C0 * h2 * h2 * det * val;
+    else
+      val = C0 * h2 * det * val;
+
+    est_el += val;
+
+    // ===== jump residuals 
+    if (C1 && (dim > 1)) {
+      int neighbours = Global::getGeo(NEIGH, dim);
+      int face;
+      Quadrature *surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
+      int numPointsSurface = surfaceQuad->getNumPoints();
+      Vector<WorldVector<double> > jump(numPointsSurface);
+
+
+      for (face = 0; face < neighbours; face++) {  
+	neigh = const_cast<Element*>(elInfo->getNeighbour(face));
+	if (neigh && neigh->getMark()) {      
+	  ElInfo *neighInfo = mesh->createNewElInfo();
+	  WorldVector<int> faceIndEl, faceIndNeigh;
+	  int oppV = elInfo->getOppVertex(face);
+	  DimVec<WorldVector<double> > LambdaNeigh(dim, NO_INIT);
+	  double detNeigh;
+	  DimVec<double> lambda(dim, NO_INIT);
+	      
+	  el->sortFaceIndices(face, &faceIndEl);
+	  neigh->sortFaceIndices(oppV, &faceIndNeigh);
+	    
+	  neighInfo->setElement(const_cast<Element*>(neigh));
+	  neighInfo->setFillFlag(Mesh::FILL_COORDS);
+	      
+	  int dow = Global::getGeo(WORLD);
+		
+	  int j, i1, i2;
+
+	  for (int i = 0; i < dow; i++)
+	    neighInfo->getCoord(oppV)[i] = elInfo->getOppCoord(face)[i];
+		
+	  // periodic leaf data ?
+	  ElementData *ldp = el->getElementData()->getElementData(PERIODIC);
+
+	  bool periodicCoords = false;
+
+	  if (ldp) {
+	    ::std::list<LeafDataPeriodic::PeriodicInfo>::iterator it;
+
+	    for (it = dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().begin();
+		 it != dynamic_cast<LeafDataPeriodic*>(ldp)->getInfoList().end();
+		 ++it) {
+	      if (it->elementSide == face) {
+		for (int i = 0; i < dim; i++) {
+		  i1 = faceIndEl[i];
+		  i2 = faceIndNeigh[i];
+
+		  for (j = 0; j < dim; j++) {
+		    if (i1 == el->getVertexOfPosition(INDEX_OF_DIM(dim - 1, 
+								   dim),
+						      face,
+						      j)) {
+		      break;
+		    }
+		  }
+
+		  TEST_EXIT_DBG(j != dim)("vertex i1 not on face ???\n");
+		      
+		  neighInfo->getCoord(i2) = (*(it->periodicCoords))[j];
+		}
+		periodicCoords = true;
+		break;
+	      }
+	    }
+	  }
+      
+	  if (!periodicCoords) {
+	    for (int i = 0; i < dim; i++) {
+	      i1 = faceIndEl[i];
+	      i2 = faceIndNeigh[i];
+	      for (j = 0; j < dow; j++)
+		neighInfo->getCoord(i2)[j] = elInfo->getCoord(i1)[j];
+	    }
+	  }
+	      
+	  Parametric *parametric = mesh->getParametric();
+	  if (parametric) {
+	    neighInfo = parametric->addParametricInfo(neighInfo);
+	  }
+	      
+	  detNeigh = abs(neighInfo->calcGrdLambda(LambdaNeigh));
+	      
+	  Vector<WorldVector<double> > grdUhEl(numPointsSurface);
+	  Vector<WorldVector<double> > grdUhNeigh(numPointsSurface);
+
+	  for (iq = 0; iq < numPointsSurface; iq++) {
+	    jump[iq].set(0.0);
+	  }
+	     
+
+	  for (int system = 0; system < numSystems; system++) {	
+	    if (matrix[system] == NULL) 
+	      continue;
+	      
+	    uh[system]->getLocalVector(el, uhEl[system]);
+		
+	    double* uhNeigh;
+	    uhNeigh = new double[basFcts[system]->getNumber()];		
+	    uh[system]->getLocalVector(neigh, uhNeigh);
+
+			
+	    for (iq = 0; iq < numPointsSurface; iq++) {
+	      lambda[face] = 0.0;
+	      for (int i = 0; i < dim; i++) {
+		lambda[faceIndEl[i]] = surfaceQuad->getLambda(iq, i);
+	      }
+		  
+	      basFcts[system]->evalGrdUh(lambda, 
+					 Lambda, 
+					 uhEl[system], 
+					 &grdUhEl[iq]);
+		  
+	      lambda[oppV] = 0.0;
+	      for (int i = 0; i < dim; i++) {
+		lambda[faceIndNeigh[i]] = surfaceQuad->getLambda(iq, i);
+	      }
+		  
+	      basFcts[system]->evalGrdUh(lambda, 
+					 LambdaNeigh, 
+					 uhNeigh, 
+					 &grdUhNeigh[iq]);
+		  
+	      grdUhEl[iq] -= grdUhNeigh[iq];
+	    }				
+
+	    delete [] uhNeigh;
+
+	    ::std::vector<double*>::iterator fac;
+
+	    for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
+		   fac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin(); 
+		 it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd(); 
+		 ++it, ++fac) {
+	      Vector<WorldVector<double> > localJump(numPointsSurface);
+	      for (iq = 0; iq < numPointsSurface; iq++) {
+		localJump[iq].set(0.0);
+	      }
+		  
+	      (*it)->weakEvalSecondOrder(numPointsSurface,
+					 grdUhEl.getValArray(),
+					 localJump.getValArray());
+	      double factor = *fac ? **fac : 1.0;
+	      if (factor != 1.0) {
+		for (int i = 0; i < numPointsSurface; i++) {
+		  localJump[i] *= factor;
+		}
+	      }
+		  
+	      for (int i = 0; i < numPointsSurface; i++) {
+		jump[i] += localJump[i];
+	      }
+	    }				     
+	  }
+	      
+  
+	  for (val = iq = 0; iq < numPointsSurface; iq++) {
+	    val += surfaceQuad->getWeight(iq) * (jump[iq] * jump[iq]);
+	  }
+	      
+	  double d = 0.5 * (det + detNeigh);
+
+	  if (norm == NO_NORM || norm == L2_NORM)
+	    val *= C1 * h2_from_det(d, dim) * d;
+	  else
+	    val *= C1 * d;
+	      
+	  if (parametric) {
+	    neighInfo = parametric->removeParametricInfo(neighInfo);
+	  }
+
+	  DELETE neighInfo;
+
+	  neigh->setEstimation(neigh->getEstimation(row) + val, row);
+
+	  est_el += val;
+
+	} 
+      } 
+       
+      val = fh[::std::max(row, 0)]->
+	getBoundaryManager()->
+	boundResidual(elInfo, matrix[::std::max(row, 0)], uh[::std::max(row, 0)]);
+      if (norm == NO_NORM || norm == L2_NORM)
+	val *= C1 * h2;
+      else
+	val *= C1;
+	
+      est_el += val;
+    } 
+  
+
+    el->setEstimation(est_el, row);
+
+    est_sum += est_el;
+    est_max = max(est_max, est_el);
+
+    elInfo->getElement()->setMark(0);   
+  }
+
+  void r(const ElInfo              *elInfo,
+	 int                        numPoints,
+	 const double              *uhIq,
+	 const WorldVector<double> *grdUhIq,
+	 const WorldMatrix<double> *D2UhIq,
+	 const double              *uhOldIq,
+	 const WorldVector<double> *grdUhOldIq,
+	 const WorldMatrix<double> *D2UhOldIq,
+	 DOFMatrix *A, 
+	 DOFVector<double> *fh,
+	 Quadrature *quad,
+	 double *result)
+  {
+    ::std::vector<Operator*>::iterator it;
+    ::std::vector<double*>::iterator fac;
+    double factor;
+
+    // lhs
+    for (it = const_cast<DOFMatrix*>(A)->getOperatorsBegin(),
+	   fac = const_cast<DOFMatrix*>(A)->getOperatorEstFactorBegin(); 
+	 it != const_cast<DOFMatrix*>(A)->getOperatorsEnd(); 
+	 ++it, ++fac) {
+      
+      factor = *fac ? **fac : 1.0;
+      if (factor) {
+	if (D2UhIq) {
+	  (*it)->evalSecondOrder(numPoints, uhIq, grdUhIq, D2UhIq, result, -factor);
+	}
+
+	if (grdUhIq) {
+	  (*it)->evalFirstOrderGrdPsi(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+	  (*it)->evalFirstOrderGrdPhi(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+	}
+	
+	if (uhIq) {
+	  (*it)->evalZeroOrder(numPoints, uhIq, grdUhIq, D2UhIq, result, factor);
+	}
+      }
+    }
+    
+    // rhs
+    for (it = const_cast<DOFVector<double>*>(fh)->getOperatorsBegin(),
+	 fac = const_cast<DOFVector<double>*>(fh)->getOperatorEstFactorBegin(); 
+	 it != const_cast<DOFVector<double>*>(fh)->getOperatorsEnd(); 
+	 ++it, ++fac) {
+
+      factor = *fac ? **fac : 1.0;
+      if (factor) {
+	if ((*it)->getUhOld()) {
+	  if (D2UhOldIq) {
+	    (*it)->evalSecondOrder(numPoints, 
+				   uhOldIq, grdUhOldIq, D2UhOldIq, 
+				   result, factor);
+	  }
+	  if (grdUhOldIq) {
+	    (*it)->evalFirstOrderGrdPsi(numPoints,