Lagrange.cc

#include "Mesh.h"
#include "Element.h"
#include "Lagrange.h"
#include "DOFAdmin.h"
#include "RCNeighbourList.h"
#include "DOFVector.h"
#include "Traverse.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "Parametric.h"

#include <stdio.h>
#include <algorithm>
#include <list>

namespace AMDiS {

  std::vector<DimVec<double>* > Lagrange::baryDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  DimVec<int>* Lagrange::ndofDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  int Lagrange::nBasFctsDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  std::vector<BasFctType*> Lagrange::phiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  std::vector<GrdBasFctType*> Lagrange::grdPhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  std::vector<D2BasFctType*> Lagrange::D2PhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  std::list<Lagrange*> Lagrange::allBasFcts;


  Lagrange::Lagrange(int dim, int degree)
    : BasisFunction(std::string("Lagrange"), dim, degree)
  {
    // set name
    char dummy[3];
    sprintf(dummy, "%d", dim);
    name.append(dummy);
    name.append(" ");
    sprintf(dummy, "%d", degree);  
    name.append(dummy);

    // set nDOF
    setNDOF();

    // set barycentric coordinates
    setBary();

    // set function pointer
    setFunctionPointer();
  }

  Lagrange::~Lagrange()
  {
    for (int i = 0; i < static_cast<int>(bary->size()); i++)
      if ((*bary)[i]) {
	delete (*bary)[i];
	(*bary)[i] = NULL;
      }
  }

  Lagrange* Lagrange::getLagrange(int dim, int degree)
  {
    std::list<Lagrange*>::iterator it;
    for (it = allBasFcts.begin(); it != allBasFcts.end(); it++)
      if (((*it)->dim == dim) && ((*it)->degree == degree))
	return (*it);

    Lagrange* newLagrange = new Lagrange(dim, degree);
    allBasFcts.push_back(newLagrange);
    return newLagrange;
  }

  void Lagrange::clear()
  {    
    for (std::list<Lagrange*>::iterator it = allBasFcts.begin(); 
	 it != allBasFcts.end(); it++)
      if (*it) {
	delete *it;
	*it = NULL;
      }     
  }

  void Lagrange::setFunctionPointer()
  {
    if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) {
      // for all positions
      for (int i = 0; i < dim+1; i++) {
	// no vertex dofs for degree 0 ?
	if (degree == 0 && i != dim) 
	  continue;
	// for all vertices/edges/...
	for (int j = 0; j < Global::getGeo(INDEX_OF_DIM(i, dim),dim); j++) {
	  // for all dofs at this position
	  for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
	    // basis function
	    phiDimDegree[dim][degree].push_back(new Phi(this, 
							INDEX_OF_DIM(i ,dim),j, k));
	    // gradients
	    grdPhiDimDegree[dim][degree].push_back(new GrdPhi(this, 
							      INDEX_OF_DIM(i, dim),
							      j, k));
	    // D2-Matrices
	    D2PhiDimDegree[dim][degree].push_back(new D2Phi(this, 
							    INDEX_OF_DIM(i, dim),
							    j, k));
	  }
	}
      }
    }
    phi = &phiDimDegree[dim][degree];
    grdPhi = &grdPhiDimDegree[dim][degree];
    d2Phi = &D2PhiDimDegree[dim][degree];

    switch (degree) {
    case 0:
      refineInter_fct = refineInter0;
      coarseRestr_fct = coarseRestr0;
      coarseInter_fct = coarseInter0;  // not yet implemented
      break;
    case 1:
      refineInter_fct = refineInter1;
      coarseRestr_fct = coarseRestr1;
      coarseInter_fct = NULL;  // not yet implemented
      break;
    case 2:
      switch (dim) {
      case 1:
	refineInter_fct = refineInter2_1d;
	coarseRestr_fct = NULL; // not yet implemented
	coarseInter_fct = coarseInter2_1d;
	break;
      case 2: 
	refineInter_fct = refineInter2_2d;
	coarseRestr_fct = coarseRestr2_2d;
	coarseInter_fct = coarseInter2_2d;
	break;
      case 3:
	refineInter_fct = refineInter2_3d;
	coarseRestr_fct = coarseRestr2_3d;
	coarseInter_fct = coarseInter2_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 3:
      switch (dim) {
      case 1:
	refineInter_fct = refineInter3_1d;
	coarseRestr_fct = coarseRestr3_1d;
	coarseInter_fct = coarseInter3_1d;
	break;
      case 2: 
	refineInter_fct = refineInter3_2d;
	coarseRestr_fct = coarseRestr3_2d;
	coarseInter_fct = coarseInter3_2d;
	break;
      case 3:
	refineInter_fct = refineInter3_3d;
	coarseRestr_fct = coarseRestr3_3d;
	coarseInter_fct = coarseInter3_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 4:
      switch (dim) {
      case 1: 
	refineInter_fct = refineInter4_1d;
	coarseRestr_fct = coarseRestr4_1d;
	coarseInter_fct = coarseInter4_1d;
	break;
      case 2:
	refineInter_fct = refineInter4_2d;
	coarseRestr_fct = coarseRestr4_2d;
	coarseInter_fct = coarseInter4_2d;
	break;
      case 3:
	refineInter_fct = refineInter4_3d;
	coarseRestr_fct = coarseRestr4_3d;
	coarseInter_fct = coarseInter4_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

  void Lagrange::setNDOF()
  {
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
      ndofDimDegree[dim][degree] = new DimVec<int>(dim, DEFAULT_VALUE, 0);

      if (degree != 0)
	(*ndofDimDegree[dim][degree])[VERTEX] = 1;
      else
	(*ndofDimDegree[dim][degree])[VERTEX] = 0;      

      for (int i = 1; i < dim + 1; i++) {
	nBasFcts = getNumberOfDOFs(i, degree);
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
	for (int j = 0; j < i; j++) {
	  (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] -= 
	    Global::getGeo(INDEX_OF_DIM(j, dim), i) * 
	    (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(j, dim)];
	}
      }
      nBasFctsDimDegree[dim][degree] = nBasFcts;
    }
    nBasFcts = nBasFctsDimDegree[dim][degree];
    nDOF = ndofDimDegree[dim][degree];
  } 

  DimVec<double> *Lagrange::getCoords(int i) const
  {
    return (*bary)[i];
  }

  void Lagrange::setVertices(int dim, int degree, 
			     GeoIndex position, int positionIndex, int nodeIndex, 
			     int** vertices)
  {
    FUNCNAME("Lagrange::setVertices()");

    TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");

    int dimOfPosition = DIM_OF_INDEX(position, dim);

    *vertices = new int[dimOfPosition + 1];

    if ((degree == 4) && (dimOfPosition==1)) {
      (*vertices)[(nodeIndex != 2) ? 0 : 1] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex,
							      0);
      (*vertices)[(nodeIndex != 2) ? 1 : 0] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex, 
							      1);
    } else if ((degree==4) && (dimOfPosition==2)) {
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+dimOfPosition*nodeIndex) % (dimOfPosition + 1)] = 
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }    
    } else {
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+nodeIndex) % (dimOfPosition + 1)] = 
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }
    }
  }

  Lagrange::Phi::Phi(Lagrange* owner, 
		     GeoIndex position, 
		     int positionIndex, 
		     int nodeIndex)
    : vertices(NULL)
  {
    FUNCNAME("Lagrange::Phi::Phi()");

    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);
  
    // set function pointer
    switch (owner->getDegree()) {
    case 0:
      switch (position) {
      case CENTER:
	func = phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
      switch (position) {
      case VERTEX:
	func = phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
      switch (position) {
      case VERTEX:
	func = phi2v;
	break;
      case EDGE:
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
	func = phi2e;
	break;
      case CENTER:
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
	func = phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
      switch (position) {
      case VERTEX:
	func = phi3v;
	break;
      case EDGE:
	func = phi3e;
	break;
      case FACE:
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
	func = phi3f;
	break;
      case CENTER:
	switch (owner->getDim()) {
	case 1:
	  func = phi3e;
	  break;
	case 2:
	  func = phi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
      switch (position) {
      case VERTEX:
	func = phi4v;
	break;
      case EDGE:
	if (nodeIndex == 1)
	  func = phi4e1;
	else 
	  func = phi4e0;
	break;
      case FACE:
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
	func = phi4f;      
	break;
      case CENTER:
	switch (owner->getDim()) {
	case 1:
	  if (nodeIndex == 1)
	    func = phi4e1;
	  else 
	    func = phi4e0;
	  break;
	case 2:
	  func = phi4f;
	  break;
	case 3:
	  func = phi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

  Lagrange::Phi::~Phi() { 
    delete [] vertices; 
  }


  Lagrange::GrdPhi::GrdPhi(Lagrange* owner, 
			   GeoIndex position, 
			   int positionIndex, 
			   int nodeIndex)
    : vertices(NULL)
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
    switch (owner->getDegree()) {
    case 0:
      switch (position) {
      case CENTER:
	func = grdPhi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
      switch (position) {
      case VERTEX:
	func = grdPhi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
      switch (position) {
      case VERTEX:
	func = grdPhi2v;
	break;
      case EDGE:
	func = grdPhi2e;
	break;
      case CENTER:
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
	func = grdPhi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
      switch (position) {
      case VERTEX:
	func = grdPhi3v;
	break;
      case EDGE:
	func = grdPhi3e;
	break;
      case FACE:
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
	func = grdPhi3f;
	break;
      case CENTER:
	switch (owner->getDim()) {
	case 1:
	  func = grdPhi3e;
	  break;
	case 2:
	  func = grdPhi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
      switch (position) {
      case VERTEX:
	func = grdPhi4v;
	break;
      case EDGE:
	if (nodeIndex == 1)
	  func = grdPhi4e1;
	else 
	  func = grdPhi4e0;
	break;
      case FACE:
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
	func = grdPhi4f;      
	break;
      case CENTER:
	switch (owner->getDim()) {
	case 1:
	  if (nodeIndex == 1)
	    func = grdPhi4e1;
	  else 
	    func = grdPhi4e0;
	  break;
	case 2:
	  func = grdPhi4f;
	  break;
	case 3:
	  func = grdPhi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

  Lagrange::GrdPhi::~GrdPhi() 
  { 
    delete [] vertices; 
  }

  Lagrange::D2Phi::D2Phi(Lagrange* owner, 
			 GeoIndex position, 
			 int positionIndex, 
			 int nodeIndex)
    : vertices(NULL)
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
    switch (owner->getDegree()) {
    case 0:
      switch (position) {
      case CENTER:
	func = D2Phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
      switch (position) {
      case VERTEX:
	func = D2Phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
      switch (position) {
      case VERTEX:
	func = D2Phi2v;
	break;
      case EDGE:
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
	func = D2Phi2e;
	break;
      case CENTER:
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
	func = D2Phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
      switch (position) {
      case VERTEX:
	func = D2Phi3v;
	break;
      case EDGE:
	func = D2Phi3e;
	break;
      case FACE:
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
	func = D2Phi3f;
	break;
      case CENTER:
	switch (owner->getDim()) {
	case 1:
	  func = D2Phi3e;
	  break;
	case 2:
	  func = D2Phi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
      switch (position) {
      case VERTEX:
	func = D2Phi4v;
	break;
      case EDGE:
	if (nodeIndex == 1)
	  func = D2Phi4e1;
	else 
	  func = D2Phi4e0;
	break;
      case FACE:
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
	func = D2Phi4f;      
	break;
      case CENTER:
	switch (owner->getDim()) {
	case 1:
	  if (nodeIndex == 1)
	    func = D2Phi4e1;
	  else 
	    func = D2Phi4e0;
	  break;
	case 2:
	  func = D2Phi4f;
	  break;
	case 3:
	  func = D2Phi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

  Lagrange::D2Phi::~D2Phi() 
  { 
    delete [] vertices; 
  }

  void Lagrange::createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
			      DimVec<double>* vec)
  {
    if (vec == NULL)
      vec = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);

    if (dimIndex == numCoords - 1) {
      (*vec)[coordInd[dimIndex]] = double(rest) / degree;
      DimVec<double>* newCoords = new DimVec<double>(*vec);
      bary->push_back(newCoords);
    } else {
      for (int i = rest - 1; i >= 1; i--) {
	(*vec)[coordInd[dimIndex]] = double(i) / degree;
	createCoords(coordInd, numCoords, dimIndex + 1, rest - i, vec);
      }
    }
  }

  void Lagrange::setBary()
  {
    bary = &baryDimDegree[dim][degree];

    if (static_cast<int>(bary->size()) == 0) {
      for (int i = 0; i <= dim; i++) { // for all positions
	int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
	for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
	  int *coordInd = new int[i + 1];      // indices of relevant coords
	  for (int k = 0; k < i + 1; k++) { 
	    coordInd[k] = Global::getReferenceElement(dim)->
	      getVertexOfPosition(INDEX_OF_DIM(i, dim), j, k);
	  }
	  createCoords(coordInd, i + 1, 0, degree);
	  delete [] coordInd;
	  if (static_cast<int>(bary->size()) == nBasFcts) 
	    return;
	}
      }
    }
  }

  int Lagrange::getNumberOfDOFs(int dim, int degree)
  {
    int result = 0;
    for (int i = 0; i <= degree; i++)
      result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));

    return result;
  }


  const DegreeOfFreedom *Lagrange::getDOFIndices(const Element *el, 
						 const DOFAdmin& admin,
						 DegreeOfFreedom *idof) const
  {
    WARNING("depricated: use getLocalIndices() instead\n");
    return getLocalIndices(el, &admin, idof);
  }

  int* Lagrange::orderOfPositionIndices(const Element* el, 
					GeoIndex position, 
					int positionIndex) const
  {
    static int sortedVertex = 0;
    static int sortedEdgeDeg2 = 0;
    static int sortedEdgeDeg3[2][2] = {{0, 1}, {1, 0}};
    static int sortedEdgeDeg4[2][3] = {{0, 1, 2}, {2, 1, 0}};
    static int sortedFaceDeg3 = 0;
    static int sortedFaceDeg4[7][3] = {{0,0,0}, {0,2,1}, {1,0,2}, {0,1,2},
				       {2,1,0}, {1,2,0}, {2,0,1}};
    static int sortedCenterDeg4 = 0;

    int dimOfPosition = DIM_OF_INDEX(position, dim);

    // vertex
    if (dimOfPosition == 0)
      return &sortedVertex;   

    // edge
    if ((dimOfPosition == 1) && (degree == 2))
      return &sortedEdgeDeg2;    
      
    int vertex[3];
    int** dof = const_cast<int**>(el->getDOF());
    int verticesOfPosition = dimOfPosition + 1;

    for (int i = 0; i < verticesOfPosition; i++)
      vertex[i] = Global::getReferenceElement(dim)->
	getVertexOfPosition(position, positionIndex, i);   
    
    if (dimOfPosition == 1) {
      if (degree == 3) {
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
	  return sortedEdgeDeg3[0];
	else
	  return sortedEdgeDeg3[1];	
      } else { // degree == 4
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
	  return sortedEdgeDeg4[0];
	else
	  return sortedEdgeDeg4[1];	
      }
    }

    // face
    if (dimOfPosition == 2) {
      if (degree == 3) {
	return &sortedFaceDeg3;
      } else {  // degree == 4!
	int no = 0;
	const Element *refElem = Global::getReferenceElement(dim);

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 0)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 1)][0]) 
	  no++;

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 1)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 2)][0])
	  no += 2;

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 2)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 0)][0])
	  no += 4;

	return sortedFaceDeg4[no];
      }
    }

    // center
    if ((dimOfPosition == 3) && (degree == 4))
      return &sortedCenterDeg4;    

    ERROR_EXIT("should not be reached\n");
    return NULL;
  }

  void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
  {
    el_info->testFlag(Mesh::FILL_BOUND);

    // boundaries
    int index = 0;
    int offset = 0;
    BoundaryType boundaryType;
    for (int i = dim - 1; i > 0; i--)
      offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);

    for (int i = 0; i < dim; i++) {
      int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
      for (int j = offset; j < jto; j++) {
	boundaryType = el_info->getBoundary(j);
	int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
	for (int k = 0; k < kto; k++) {
	  bound[index++] = boundaryType;
	}
      }
      offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
    }

    // interior nodes in the center
    for (int i = 0; i < (*nDOF)[CENTER]; i++) {
      bound[index++] = INTERIOR;
    }

    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
  }


  const double* Lagrange::interpol(const ElInfo *el_info, 
				   int no, const int *b_no, 
				   AbstractFunction<double, WorldVector<double> > *f, 
				   double *vec)
  {
    FUNCNAME("Lagrange::interpol()");

    static double* localVec = NULL;
    static int localVecSize = 0;

    double *rvec = NULL;

    if (vec) {
      rvec = vec;
    } else {
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new double[nBasFcts]; 
      }
      if (!localVec)
	localVec = new double[nBasFcts]; 

      localVecSize = nBasFcts;
      rvec = localVec;
    } 

    WorldVector<double> x;

    el_info->testFlag(Mesh::FILL_COORDS);
 
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	rvec[0] = 0.0;
	return(const_cast<const double *>( rvec));
      }
	
      for (int i = 0; i < no; i++) {
	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
	} else {
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
	  rvec[i] = (*f)(x);
	}
      }
    } else {
      int vertices = Global::getGeo(VERTEX, dim);
      for (int i = 0; i < vertices; i++)
	rvec[i] = (*f)(el_info->getCoord(i));
      for (int i = vertices; i < nBasFcts; i++) {
	el_info->coordToWorld(*(*bary)[i], x);
	rvec[i] = (*f)(x);
      }
    }  
    
    return(const_cast<const double *>( rvec));
  }

  const WorldVector<double>* 
  Lagrange::interpol(const ElInfo *el_info, int no, 
		     const int *b_no,
		     AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
		     WorldVector<double> *vec)
  {
    FUNCNAME("*Lagrange::interpol_d()");

    static WorldVector<double> *inter;
    WorldVector<double> *rvec = 
      vec ? vec : (inter ? inter : inter = new WorldVector<double>[getNumber()]);
    WorldVector<double> x;

    el_info->testFlag(Mesh::FILL_COORDS);
  
    int vertices = Global::getGeo(VERTEX, dim);

    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	for (int i = 0; i < dow; i++)
	  rvec[0][i] = 0.0;
	return(const_cast<const WorldVector<double> *>( rvec));
      }
      for (int i = 0; i < no; i++) {
	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
	} else {
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
	  rvec[i] = (*f)(x);
	}
      }
    } else {
      for (int i = 0; i < vertices; i++)
	rvec[i] = (*f)(el_info->getCoord(i));
      for (int i = vertices; i < nBasFcts; i++) {
	el_info->coordToWorld(*(*bary)[i], x);
	rvec[i] = (*f)(x);
      }
    }  
    
    return(const_cast<const WorldVector<double> *>( rvec));
  }


  const DegreeOfFreedom *Lagrange::getLocalIndices(const Element* el,
						   const DOFAdmin *admin,
						   DegreeOfFreedom *indices) const
  {
    static DegreeOfFreedom *localVec = NULL;
    static int localVecSize = 0;

    const DegreeOfFreedom **dof = el->getDOF();

    int nrDOFs, n0, node0, num = 0;
    GeoIndex posIndex;
    DegreeOfFreedom* result;

    if (indices) {
      result = indices;
    } else {
      if (localVec && nBasFcts > localVecSize) {
	delete [] localVec;
	localVec = new DegreeOfFreedom[nBasFcts];
      }
      if (!localVec)
	localVec = new DegreeOfFreedom[nBasFcts];

      localVecSize = nBasFcts;
      result = localVec;
    }

    for (int pos = 0, j = 0; pos <= dim; pos++) {
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

	for (int i = 0; i < num; node0++, i++) {
	  const int *indi = orderOfPositionIndices(el, posIndex, i);

	  for (int k = 0; k < nrDOFs; k++)
	    result[j++] = dof[node0][n0 + indi[k]];  
	}
      }
    }

    return result;
  }

  void Lagrange::getLocalIndicesVec(const Element* el,
				    const DOFAdmin *admin,
				    Vector<DegreeOfFreedom> *indices) const
  {
    if (indices->getSize() < nBasFcts)
      indices->resize(nBasFcts);

    getLocalIndices(el, admin, &((*indices)[0]));
  }

  void Lagrange::getLocalDofPtrVec(const Element *el, 
				   const DOFAdmin *admin,
				   std::vector<const DegreeOfFreedom*>& vec) const
  {
    if (static_cast<int>(vec.size()) < nBasFcts)
      vec.resize(nBasFcts);

    const DegreeOfFreedom **dof = el->getDOF();
    GeoIndex posIndex;
    int nrDOFs, n0, node0, num = 0;

    for (int pos = 0, j = 0; pos <= dim; pos++) {
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

	for (int i = 0; i < num; node0++, i++) {
	  const int *indi = orderOfPositionIndices(el, posIndex, i);

	  for (int k = 0; k < nrDOFs; k++)
	    vec[j++] = &(dof[node0][n0 + indi[k]]);
	}
      }
    }
  }

  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
    ERROR_EXIT("not yet\n");
  }

  void Lagrange::l2ScpFctBas(Quadrature *quad,
			     AbstractFunction<double, WorldVector<double> >* f,
			     DOFVector<double>* fh)
  {
    FUNCNAME("Lagrange::l2ScpFctBas()");

    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
    TEST_EXIT_DBG(fh->getFESpace())
      ("no fe_space in DOF_REAL_VEC %s\n", fh->getName().c_str());
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)
      ("wrong basis fcts for fh\n");
    TEST_EXIT_DBG(!fh->getFESpace()->getMesh()->getParametric())
      ("Not yet implemented!");

    if (!quad)
      quad = Quadrature::provideQuadrature(dim, 2 * degree - 2);

    const FastQuadrature *quad_fast = 
      FastQuadrature::provideFastQuadrature(this, *quad, INIT_PHI);
    double *wdetf_qp = new double[quad->getNumPoints()];
    int nPoints = quad->getNumPoints();
    DOFAdmin *admin = fh->getFESpace()->getAdmin();
    WorldVector<double> x;

    TraverseStack stack;
    ElInfo *el_info = stack.traverseFirst(fh->getFESpace()->getMesh(), -1, 
					  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (el_info) {
      const DegreeOfFreedom  *dof = getLocalIndices(el_info->getElement(), admin, NULL);
      double det = el_info->getDet();

      for (int iq = 0; iq < nPoints; iq++) {
	el_info->coordToWorld(quad->getLambda(iq), x);
	wdetf_qp[iq] = quad->getWeight(iq) * det * ((*f)(x));
      }
      
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
	for (int iq = 0; iq < nPoints; iq++)
	  val += quad_fast->getPhi(iq, j) * wdetf_qp[iq];

	(*fh)[dof[j]] += val;
      }

      el_info = stack.traverseNext(el_info);
    }

    delete [] wdetf_qp;
  }


  // ===== refineInter functions =====

  void  Lagrange::refineInter0(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter0()");

    if (n < 1) 
      return;

    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(CENTER);
    Element *el = list->getElement(0);
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);
    // Parent center
    DegreeOfFreedom dof0 = el->getDOF(node, n0);           
    // Newest vertex is center 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0);  

    (*drv)[dof_new] = (*drv)[dof0];
    // Newest vertex is center 
    dof_new = el->getChild(1)->getDOF(node, n0);  
    (*drv)[dof_new] = (*drv)[dof0];
  }

  void Lagrange::refineInter1(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
			      BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter1_1d()");

    if (n < 1) 
      return;

    int dim = drv->getFESpace()->getMesh()->getDim();
    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
    Element *el = list->getElement(0);
    // 1st endpoint of refinement edge
    DegreeOfFreedom dof0 = el->getDOF(0, n0);
    // 2nd endpoint of refinement edge          
    DegreeOfFreedom dof1 = el->getDOF(1, n0);
    // newest vertex is DIM 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0); 
    (*drv)[dof_new] = 0.5 * ((*drv)[dof0] + (*drv)[dof1]);
  }

  void Lagrange::refineInter2_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
				 BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter2_1d()");

    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
  
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);

    /****************************************************************************/
    /*  newest vertex of child[0] and child[1]                                  */
    /****************************************************************************/

    // newest vertex is DIM
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node + 1, n0);
    (*drv)[cdof] = (*drv)[pdof[2]];

    node = drv->getFESpace()->getMesh()->getNode(CENTER);        
    n0 = admin->getNumberOfPreDOFs(CENTER);

    /****************************************************************************/
    /*  midpoint of edge on child[0] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(0)->getDOF(node, n0); 
    (*drv)[cdof] = 
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node, n0); 
    (*drv)[cdof] = 
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
  }

  void Lagrange::refineInter2_2d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter2_2d()");

    if (n < 1) 
      return;

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
  
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);

    /****************************************************************************/
    /*  newest vertex of child[0] and child[1]                                  */
    /****************************************************************************/

    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
    (*drv)[cdof] = (*drv)[pdof[5]];

    node = drv->getFESpace()->getMesh()->getNode(EDGE);        
    n0 = admin->getNumberOfPreDOFs(EDGE);

    /****************************************************************************/
    /*  midpoint of edge on child[0] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(0)->getDOF(node, n0); 
    (*drv)[cdof] = 
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];

    /****************************************************************************/
    /* node in the common edge of child[0] and child[1]                         */
    /****************************************************************************/

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]]  + 0.75 * (*drv)[pdof[5]];

    if (n > 1) {
      /****************************************************************************/
      /*  adjust the value at the midpoint of the common edge of neigh's children */
      /****************************************************************************/
      el = list->getElement(1);
      pdof = basFct->getLocalIndices(el, admin, NULL);
      
      cdof = el->getChild(0)->getDOF(node+1, n0); 
      (*drv)[cdof] = 
	-0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
	+ 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
    }
  }

  void Lagrange::refineInter2_3d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter2_3d()");

    if (n < 1) 
      return;

    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();

    DegreeOfFreedom pdof[10];
    basFct->getLocalIndices(el, admin, pdof);

    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);

    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/

    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);

    (*drv)[cdof[3]] = ((*drv)[pdof[4]]);
    (*drv)[cdof[6]] = 
      (0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
    (*drv)[cdof[8]] = 
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
    (*drv)[cdof[9]] = 
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/
  
    DegreeOfFreedom cdofi = el->getChild(1)->getDOF(node0 + 2, n0);
    (*drv)[cdofi] = 
      (-0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
    for (int i = 1; i < n; i++) {
      el = list->getElement(i);
      basFct->getLocalIndices(el, admin, pdof);
      
      int lr_set = 0;
      if (list->getNeighbourElement(i, 0) && list->getNeighbourNr(i, 0) < i)
	lr_set = 1;
      
      if (list->getNeighbourElement(i, 1) && list->getNeighbourNr(i, 1) < i)
	lr_set += 2;
      
      TEST_EXIT(lr_set)("no values set on both neighbours\n");
      
      /****************************************************************************/
      /*  values on child[0]                                                      */
      /****************************************************************************/
      
      switch (lr_set) {
      case 1:
	cdofi = el->getChild(0)->getDOF(node0+4, n0);      
	(*drv)[cdofi] = 
	  (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
	   + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
	break;
      case 2:
	cdofi = el->getChild(0)->getDOF(node0+5, n0);      
	(*drv)[cdofi] = 
	  (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
	   + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
	
      }
    }
  }

  void Lagrange::refineInter3_1d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter3_1d()");

    if (n < 1) 
      return;

    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();

    DegreeOfFreedom *pdof = new DegreeOfFreedom[basFct->getNumber()];
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);

    (*drv)[cdof[2]] = 
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
    (*drv)[cdof[3]] = 
      (0.3125 * ((*drv)[pdof[0]] - (*drv)[pdof[8]]) + 0.0625 * (*drv)[pdof[1]]
       + 0.9375 * (*drv)[pdof[7]]);
    (*drv)[cdof[4]] = (*drv)[pdof[7]];
    (*drv)[cdof[5]] = (*drv)[pdof[9]];
    (*drv)[cdof[6]] =  
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
    (*drv)[cdof[9]] = 
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.125 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

    cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
    (*drv)[cdof[9]] =  
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);

    if (n <= 1) {
      delete [] pdof;
      return;
    }
    /****************************************************************************/
    /*  adjust the value on the neihgbour                                       */
    /****************************************************************************/

    el = list->getElement(1);
    basFct->getLocalIndices(el, admin, pdof);

    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/
    cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
  
    (*drv)[cdof[5]] =  (*drv)[pdof[9]];
    (*drv)[cdof[6]] = 
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
    (*drv)[cdof[9]] = 
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.12500 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);

    (*drv)[dof9] =  
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) +  0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75  *(*drv)[pdof[9]]);

    delete [] pdof;
  }

  void Lagrange::refineInter3_2d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, 
				 BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter3_2d()");

    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();

    DegreeOfFreedom pdof[10];
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);

    (*drv)[cdof[2]] =
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
    (*drv)[cdof[3]] = 
      (0.3125 * ((*drv)[pdof[0]] - (*drv)[pdof[8]]) + 0.0625 * (*drv)[pdof[1]]
       + 0.9375 * (*drv)[pdof[7]]);
    (*drv)[cdof[4]] = (*drv)[pdof[7]];
    (*drv)[cdof[5]] = (*drv)[pdof[9]];
    (*drv)[cdof[6]] =  
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
    (*drv)[cdof[9]] = 
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.125 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);

    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/

    cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);

    (*drv)[cdof[5]] = (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));