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()
  {
  }

  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::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)
  {
    TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");
    int dimOfPosition = DIM_OF_INDEX(position, dim);

    *vertices = GET_MEMORY(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 = GET_MEMORY(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);
	  FREE_MEMORY(coordInd, int, i + 1);
	  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;
  }

  const BoundaryType* Lagrange::getBound(const ElInfo* el_info, 
					 BoundaryType* bound) const
  {
    static BoundaryType *bound_vec = NULL;
    static int bound_vec_size = 0;
    BoundaryType *result;

    if (bound) {
      result = bound;
    } else {
      // realloc memory if neccessary
      if (bound_vec_size < nBasFcts) {
	if (bound_vec) 
	  FREE_MEMORY(bound_vec, BoundaryType, bound_vec_size);
	bound_vec = GET_MEMORY(BoundaryType, nBasFcts);
	bound_vec_size = nBasFcts;
      }

      result = bound_vec;
    }

    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++) {
	  result[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++) {
      result[index++] = INTERIOR;
    }

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

    return result;
  }


  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 *inter = NULL;
    static int inter_size = 0;

    double *rvec =  NULL;

    if (vec) {
      rvec = vec;
    } else {
      if(inter) FREE_MEMORY(inter, double, nBasFcts);
      inter = GET_MEMORY(double, nBasFcts);
      inter_size = nBasFcts;
      rvec = inter;
    } 

    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();
    const int *indi;
    int nrDOFs, n0, node0, num = 0;
    GeoIndex posIndex;

    DegreeOfFreedom* result;
    
    if (indices) {
      result = indices;
    } else {
      if (localVec) 
	FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
      localVec = GET_MEMORY(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++) {
	  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::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
    ERROR_EXIT("not yet\n");
  }

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

    TraverseStack stack;
    double *wdetf_qp = NULL;
    double val, det;
    WorldVector<double> x;
    int iq, j;
    const DegreeOfFreedom  *dof;
    const Parametric *parametric;

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