ElInfo3d.cc

#include "ElInfo3d.h"
#include "BasisFunction.h"
#include "Element.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "FiniteElemSpace.h"
#include "Flag.h"
#include "MacroElement.h"
#include "Mesh.h"
#include "Global.h"
#include "FixVec.h"
#include "DOFVector.h"

namespace AMDiS {

  void ElInfo3d::fillMacroInfo(const MacroElement * mel)
  {
    FUNCNAME("ElInfo3d::fillMacroInfo()");
    Element *nb;
    MacroElement *mnb;
    Flag fill_opp_coords;

    macroElement  = const_cast<MacroElement*>( mel);
    element  = const_cast<Element*>( mel->getElement());
    parent = NULL;
    level = 0;
    elType = const_cast<MacroElement*>(mel)->getElType();

    int vertices = mesh->getGeo(VERTEX);

    if (fillFlag.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
      for (int i = 0; i < vertices; i++)
	coord[i] = mel->coord[i];
    }

    int neighbours = mesh->getGeo(NEIGH);

    if (fillFlag.isSet(Mesh::FILL_OPP_COORDS) || 
        fillFlag.isSet(Mesh::FILL_NEIGH)) {

      fill_opp_coords.setFlags(fillFlag & Mesh::FILL_OPP_COORDS);
      for (int i = 0; i < neighbours; i++) {
	if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
	  neighbour[i] = const_cast<Element*>( mel->getNeighbour(i)->getElement());
	  nb = const_cast<Element*>( neighbour[i]);
	  int k;
	  k = oppVertex[i] = mel->getOppVertex(i);

	  if (nb->getChild(0) && (k < 2)) {   /*make nb nearest element.*/
	    if (k == 1) {
	      neighbour[i]      = const_cast<Element*>( nb->getChild(0));
	      nb = const_cast<Element*>( neighbour[i]);
	    } else {
	      neighbour[i]      = const_cast<Element*>( nb->getChild(1));
	      nb = const_cast<Element*>( neighbour[i]);
	    }
	    k = oppVertex[i] = 3;
	    if (fill_opp_coords.isAnySet()) {
	      /* always edge between vertices 0 and 1 is bisected! */
	      if (mnb->getElement()->isNewCoordSet())
		oppCoord[i] = *(mnb->getElement()->getNewCoord());
	      else
		oppCoord[i] = (mnb->coord[0] + mnb->coord[1]) * 0.5;
	    }
	  } else {
	    if  (fill_opp_coords.isAnySet()) {
	      oppCoord[i] = mnb->coord[k];
	    }
	  }
	} else {
	  neighbour[i] = NULL;
	}
      }
    }

    if (fillFlag.isSet(Mesh::FILL_BOUND)) {
      for (int i = 0; i < element->getGeo(BOUNDARY); i++)
	boundary[i] = mel->getBoundary(i);     

      for (int i = 0; i < element->getGeo(PROJECTION); i++)
	projection[i] = mel->getProjection(i);
    }

    if (fillFlag.isSet(Mesh::FILL_ORIENTATION)) {
      WorldVector<WorldVector<double> > a;
      double s;

      for (int i = 0; i < 3; i++) {
	a[i] = mel->coord[i + 1];
	a[i] -= mel->coord[0];
      }

      s = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) * a[2][0]
	+ (a[0][2] * a[1][0] - a[0][0] * a[1][2]) * a[2][1]
	+ (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * a[2][2];

      if (s >= 0)
	orientation = 1;
      else
	orientation = -1;
    }
  }

  double ElInfo3d::calcGrdLambda(DimVec<WorldVector<double> >& grd_lam)
  {
    FUNCNAME("ElInfo3d::calcGrdLambda()");

    TEST_EXIT_DBG(dimOfWorld == 3)
      ("dim != dim_of_world ! use parametric elements!\n");

    WorldVector<double> *e1 = &tmpWorldVecs[0];
    WorldVector<double> *e2 = &tmpWorldVecs[1];
    WorldVector<double> *e3 = &tmpWorldVecs[2];
    WorldVector<double> *v0 = &tmpWorldVecs[3];

    double det, adet;
    double a11, a12, a13, a21, a22, a23, a31, a32, a33;

    testFlag(Mesh::FILL_COORDS);

    *v0 = coord[0];
    for (int i = 0; i < 3; i++) {
      (*e1)[i] = coord[1][i] - (*v0)[i];
      (*e2)[i] = coord[2][i] - (*v0)[i];
      (*e3)[i] = coord[3][i] - (*v0)[i];
    }

    det = (*e1)[0] * ((*e2)[1] * (*e3)[2] - (*e2)[2] * (*e3)[1])
      - (*e1)[1] * ((*e2)[0] * (*e3)[2] - (*e2)[2] * (*e3)[0])
      + (*e1)[2] * ((*e2)[0] * (*e3)[1] - (*e2)[1] * (*e3)[0]);

    adet = abs(det);

    if (adet < 1.0E-25) {
      MSG("abs(det) = %f\n",adet);
      for (int i = 0; i < 4; i++)
	for (int j = 0; j < 3; j++)
	  grd_lam[i][j] = 0.0;
    } else {
      det = 1.0 / det;
      a11 = ((*e2)[1] * (*e3)[2] - (*e2)[2] * (*e3)[1]) * det;    /* (a_ij) = A^{-T} */
      a12 = ((*e2)[2] * (*e3)[0] - (*e2)[0] * (*e3)[2]) * det;
      a13 = ((*e2)[0] * (*e3)[1] - (*e2)[1] * (*e3)[0]) * det;
      a21 = ((*e1)[2] * (*e3)[1] - (*e1)[1] * (*e3)[2]) * det;
      a22 = ((*e1)[0] * (*e3)[2] - (*e1)[2] * (*e3)[0]) * det;
      a23 = ((*e1)[1] * (*e3)[0] - (*e1)[0] * (*e3)[1]) * det;
      a31 = ((*e1)[1] * (*e2)[2] - (*e1)[2] * (*e2)[1]) * det;
      a32 = ((*e1)[2] * (*e2)[0] - (*e1)[0] * (*e2)[2]) * det;
      a33 = ((*e1)[0] * (*e2)[1] - (*e1)[1] * (*e2)[0]) * det;

      grd_lam[1][0] = a11;
      grd_lam[1][1] = a12;
      grd_lam[1][2] = a13;
      grd_lam[2][0] = a21;
      grd_lam[2][1] = a22;
      grd_lam[2][2] = a23;
      grd_lam[3][0] = a31;
      grd_lam[3][1] = a32;
      grd_lam[3][2] = a33;

      grd_lam[0][0] = -grd_lam[1][0] - grd_lam[2][0] - grd_lam[3][0];
      grd_lam[0][1] = -grd_lam[1][1] - grd_lam[2][1] - grd_lam[3][1];
      grd_lam[0][2] = -grd_lam[1][2] - grd_lam[2][2] - grd_lam[3][2];
    }

    return adet;
  }

  const int ElInfo3d::worldToCoord(const WorldVector<double>& xy,
				   DimVec<double>* lambda) const
  {
    FUNCNAME("ElInfo::worldToCoord()");

    DimVec<WorldVector<double> > edge(mesh->getDim(), NO_INIT);
    WorldVector<double> x;
    double  x0, det, det0, det1, det2;
  
    static DimVec<double> vec(mesh->getDim(), NO_INIT);

    TEST_EXIT_DBG(lambda)("lambda must not be NULL\n");

    int dim = mesh->getDim();

    TEST_EXIT_DBG(dim == dimOfWorld)("dim!=dimOfWorld not yet implemented\n");
    
    /*  wir haben das gleichungssystem zu loesen: */
    /*       ( q1x q2x q3x)  (lambda1)     (qx)      */
    /*       ( q1y q2y q3y)  (lambda2)  =  (qy)      */
    /*       ( q1z q2z q3z)  (lambda3)     (qz)      */
    /*      mit qi=pi-p3, q=xy-p3                 */

    for (int j = 0; j < dimOfWorld; j++) {
      x0 = coord[dim][j];
      x[j] = xy[j] - x0;

      for (int i = 0; i < dim; i++)
	edge[i][j] = coord[i][j] - x0;
    }

    det =  edge[0][0] * edge[1][1] * edge[2][2]
      + edge[0][1] * edge[1][2] * edge[2][0]
      + edge[0][2] * edge[1][0] * edge[2][1]
      - edge[0][2] * edge[1][1] * edge[2][0]
      - edge[0][0] * edge[1][2] * edge[2][1]
      - edge[0][1] * edge[1][0] * edge[2][2];
    det0 =       x[0] * edge[1][1] * edge[2][2]
      +       x[1] * edge[1][2] * edge[2][0]
      +       x[2] * edge[1][0] * edge[2][1]
      -       x[2] * edge[1][1] * edge[2][0]
      -       x[0] * edge[1][2] * edge[2][1]
      -       x[1] * edge[1][0] * edge[2][2];
    det1 = edge[0][0] *       x[1] * edge[2][2]
      + edge[0][1] *       x[2] * edge[2][0]
      + edge[0][2] *       x[0] * edge[2][1]
      - edge[0][2] *       x[1] * edge[2][0]
      - edge[0][0] *       x[2] * edge[2][1]
      - edge[0][1] *       x[0] * edge[2][2];
    det2 = edge[0][0] * edge[1][1] *       x[2]
      + edge[0][1] * edge[1][2] *       x[0]
      + edge[0][2] * edge[1][0] *       x[1]
      - edge[0][2] * edge[1][1] *       x[0]
      - edge[0][0] * edge[1][2] *       x[1]
      - edge[0][1] * edge[1][0] *       x[2];
  
    if (abs(det) < DBL_TOL) {
      ERROR("det = %le; abort\n", det);

      for (int i = 0; i <= dim; i++)
	(*lambda)[i] = 1.0 / dim;

      return 0;
    }

    (*lambda)[0] = det0 / det;
    (*lambda)[1] = det1 / det;
    (*lambda)[2] = det2 / det;
    (*lambda)[3] = 1.0 - (*lambda)[0] - (*lambda)[1] - (*lambda)[2];
  
    int k = -1;
    double lmin = 0.0;

    for (int i = 0; i <= dim; i++) {
      if ((*lambda)[i] < -1.E-5) {
	if ((*lambda)[i] < lmin) {
	  k = i;
	  lmin = (*lambda)[i];
	}
      }
    }

    return k;
  }


  /****************************************************************************/
  /*   update EL_INFO structure after refinement (of some neighbours)	    */
  /****************************************************************************/

  void ElInfo3d::update()
  {
    FUNCNAME("ElInfo::update()");

    int neighbours = mesh->getGeo(NEIGH);
    int vertices = mesh->getGeo(VERTEX);
  
    if (fillFlag.isSet(Mesh::FILL_NEIGH) || fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {
      Tetrahedron *nb;
      Flag fill_opp_coords = fillFlag & Mesh::FILL_OPP_COORDS;
      
      for (int ineigh = 0; ineigh < neighbours; ineigh++) {
	if ((nb = dynamic_cast<Tetrahedron*>(const_cast<Element*>(neighbour[ineigh])))) {
	  int ov = oppVertex[ineigh];
	  if (ov < 2 && nb->getFirstChild()) {
	    if (fill_opp_coords != Flag(0)) {
	      int k = -1;
	      for (int j = 0; j < vertices; j++)
		if (element->getDOF(j) == nb->getDOF(1 - ov)) 
		  k = j;
	      
	      if (k == -1) {
		for (int j = 0; j < vertices; j++)
		  if (mesh->associated(element->getDOF(j, 0), nb->getDOF(1 - ov, 0)))
		    k = j;
	      }
	      TEST_EXIT_DBG(k >= 0)("neighbour dof not found\n");
	      
	      if (nb->isNewCoordSet())
		oppCoord[ineigh] = *(nb->getNewCoord());
	      else
		for (int j = 0; j < dimOfWorld; j++)
		  oppCoord[ineigh][j] = (oppCoord[ineigh][j] + coord[k][j]) / 2;
	    }
	    neighbour[ineigh] = dynamic_cast<Tetrahedron*>(const_cast<Element*>(nb->getChild(1-ov)));
	    oppVertex[ineigh] = 3;
	  }
	}
      }
    }
  }


  double ElInfo3d::getNormal(int face, WorldVector<double> &normal)
  {
    FUNCNAME("ElInfo3d::getNormal()");

    double det = 0.0;

    WorldVector<double> *e0 = &tmpWorldVecs[0];
    WorldVector<double> *e1 = &tmpWorldVecs[1];
    WorldVector<double> *e2 = &tmpWorldVecs[2];

    if (dimOfWorld == 3) {
      int i0 = (face + 1) % 4;
      int i1 = (face + 2) % 4;
      int i2 = (face + 3) % 4;

      for (int i = 0; i < dimOfWorld; i++) {
	(*e0)[i] = coord[i1][i] - coord[i0][i];
	(*e1)[i] = coord[i2][i] - coord[i0][i];
	(*e2)[i] = coord[face][i] - coord[i0][i];
      }

      vectorProduct(*e0, *e1, normal);

      if ((*e2 * normal) < 0.0)
	for (int i = 0; i < dimOfWorld; i++)
	  normal[i] = -normal[i];

      det = norm(&normal);
      TEST_EXIT_DBG(det > 1.e-30)("det = 0 on face %d\n", face);

      normal[0] /= det;
      normal[1] /= det;
      normal[2] /= det;
    } else {
      MSG("not implemented for DIM_OF_WORLD = %d in 3d\n", dimOfWorld);
    }

    return(det);
  }




  void ElInfo3d::fillElInfo(int ichild, const ElInfo *elInfoOld)
  {
    FUNCNAME("ElInfo3d::fillElInfo()");

    int ochild = 0;             /* index of other child = 1-ichild */
    int *cv = NULL;             /* cv = child_vertex[el_type][ichild] */
    const int (*cvg)[4] = NULL;     /* cvg = child_vertex[el_type] */
    int *ce;                    /* ce = child_edge[el_type][ichild] */
    Element *nb, *nbk;
    Element *el_old = elInfoOld->element;
    Flag fillFlag_local = elInfoOld->fillFlag;
    DegreeOfFreedom *dof;
    int ov = -1;
    Mesh *mesh = elInfoOld->getMesh();

    TEST_EXIT_DBG(el_old->getChild(0))("missing child?\n"); 

    element = const_cast<Element*>( el_old->getChild(ichild));
    macroElement = elInfoOld->macroElement;
    fillFlag = fillFlag_local;
    parent  = el_old;
    level = elInfoOld->level + 1;
    iChild = ichild;
    int el_type_local = (dynamic_cast<const ElInfo3d*>(elInfoOld))->getType();
    elType = (el_type_local + 1) % 3;

    TEST_EXIT_DBG(element)("missing child %d?\n", ichild);

    if (fillFlag_local.isAnySet()) {
      cvg = Tetrahedron::childVertex[el_type_local];
      cv = const_cast<int*>(cvg[ichild]);
      ochild = 1 - ichild;
    }

    if (fillFlag_local.isSet(Mesh::FILL_COORDS) || 
	fillFlag.isSet(Mesh::FILL_DET) ||
	fillFlag.isSet(Mesh::FILL_GRD_LAMBDA)) {
      for (int i = 0; i < 3; i++) {
	for (int j = 0; j < dimOfWorld; j++)
	  coord[i][j] = elInfoOld->coord[cv[i]][j];
      }
      if (el_old->getNewCoord()) {
	coord[3] = *(el_old->getNewCoord());
      } else {
	for (int j = 0; j < dimOfWorld; j++)
	  coord[3][j] = (elInfoOld->coord[0][j] + elInfoOld->coord[1][j]) / 2;
      }
    }

    if (fillFlag_local.isSet(Mesh::FILL_NEIGH) || 
	fillFlag.isSet(Mesh::FILL_OPP_COORDS)) {

      FixVec<Element*, NEIGH> *neigh_local = &neighbour;
      const FixVec<Element*, NEIGH> *neigh_old = &elInfoOld->neighbour;

      Flag fill_opp_coords;
      fill_opp_coords.setFlags(fillFlag_local & Mesh::FILL_OPP_COORDS);
      
      /*----- nb[0] is other child --------------------------------------------*/
      
      /*    if (nb = el_old->child[ochild])   old version  */
      if (el_old->getChild(0)  &&  
	  (nb = const_cast<Element*>( el_old->getChild(ochild)))) {
	
	if (nb->getChild(0)) {         /* go down one level for direct neighbour */
	  if (fill_opp_coords.isAnySet()) {
	    if (nb->getNewCoord()) {
	      oppCoord[0]= *(nb->getNewCoord());
	    } else {
	      int k = cvg[ochild][1];
	      for (int j = 0; j < dimOfWorld; j++)
		oppCoord[0][j] = (elInfoOld->coord[ochild][j] + elInfoOld->coord[k][j]) / 2;
	    }
	  }
	  (*neigh_local)[0] = const_cast<Element*>( nb->getChild(1));
	  oppVertex[0] = 3;
	} else {
	  if (fill_opp_coords.isAnySet())
	    for (int j = 0; j < dimOfWorld; j++)
	      oppCoord[0][j] = elInfoOld->coord[ochild][j];

	  (*neigh_local)[0] = nb;
	  oppVertex[0] = 0;
	}
      } else {
	ERROR_EXIT("no other child");
	(*neigh_local)[0] = NULL;
      }


      /*----- nb[1],nb[2] are childs of old neighbours nb_old[cv[i]] ----------*/
      
      for (int i = 1; i < 3; i++) {
	if ((nb = const_cast<Element*>( (*neigh_old)[cv[i]]))) {
	  TEST_EXIT_DBG(nb->getChild(0))("nonconforming triangulation\n");
	  
	  int k;
	  for (k = 0; k < 2; k++) { /* look at both childs of old neighbour */
	    nbk = const_cast<Element*>( nb->getChild(k));
	    if (nbk->getDOF(0) == el_old->getDOF(ichild)) {
	      /* opp. vertex */
	      dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex[cv[i]])); 
	      
	      if (dof == nbk->getDOF(1)) {
		ov = 1;
		if (nbk->getChild(0)) {
		  if (fill_opp_coords.isAnySet()) {
		    if (nbk->getNewCoord()) {
		      oppCoord[i] = *(nbk->getNewCoord());
		    } else {
		      for (int j = 0; j < dimOfWorld; j++)
			oppCoord[i][j] = (elInfoOld->oppCoord[cv[i]][j] + 
					   elInfoOld->coord[ichild][j]) / 2;		      
		    }
		  }
		  (*neigh_local)[i] = nbk->getChild(0);
		  oppVertex[i] = 3;
		  break;
		}
	      } else {
		if (dof != nbk->getDOF(2)) { 
		  ov = -1; 
		  break; 
		}
		ov = 2;
	      }

	      if (fill_opp_coords.isAnySet())
		for (int j = 0; j < dimOfWorld; j++)
		  oppCoord[i][j] = elInfoOld->oppCoord[cv[i]][j];

	      (*neigh_local)[i] = nbk;
	      oppVertex[i] = ov;
	      break;
	    }
	    
	  } /* end for k */

	  // periodic ?
	  if (k == 2 || ov == -1) {
	    for (k = 0; k < 2; k++) {  /* look at both childs of old neighbour */
	      nbk = const_cast<Element*>( nb->getChild(k));
	      if (nbk->getDOF(0) == el_old->getDOF(ichild) ||
		  mesh->associated(nbk->getDOF(0, 0), el_old->getDOF(ichild, 0))) {
		/* opp. vertex */
		dof = const_cast<int*>(nb->getDOF(elInfoOld->oppVertex[cv[i]])); 
		
		if (dof == nbk->getDOF(1) || 
		    mesh->associated(dof[0], nbk->getDOF(1, 0))) {
		  ov = 1;
		  if (nbk->getChild(0)) {
		    if (fill_opp_coords.isAnySet()) {
		      if (nbk->getNewCoord()) {
			oppCoord[i] = *(nbk->getNewCoord());
		      } else {
			for (int j = 0; j < dimOfWorld; j++)
			  oppCoord[i][j] = (elInfoOld->oppCoord[cv[i]][j] + 
					    elInfoOld->coord[ichild][j]) / 2;
		      }
		    }
		    (*neigh_local)[i] = nbk->getChild(0);
		    oppVertex[i] = 3;
		    break;
		  }
		} else {
		  TEST_EXIT_DBG(dof == nbk->getDOF(2) || 
				mesh->associated(dof[0], nbk->getDOF(2, 0)))
		    ("opp_vertex not found\n");
		  ov = 2;
		}

		if (fill_opp_coords.isAnySet())
		  for (int j = 0; j < dimOfWorld; j++)
		    oppCoord[i][j] = elInfoOld->oppCoord[cv[i]][j];

		(*neigh_local)[i] = nbk;
		oppVertex[i] = ov;
		break;
	      }
	      
	    } /* end for k */

	    TEST_EXIT_DBG(k < 2)("child not found with vertex\n");
	  }
	} else {
	  (*neigh_local)[i] = NULL;
	}
      }  /* end for i */
      
      
      /*----- nb[3] is old neighbour neigh_old[ochild] ------------------------*/
      
      if (((*neigh_local)[3] = (*neigh_old)[ochild])) {
	oppVertex[3] = elInfoOld->oppVertex[ochild];

	if (fill_opp_coords.isAnySet())
	  for (int j = 0; j < dimOfWorld; j++)
	    oppCoord[3][j] = elInfoOld->oppCoord[ochild][j];
      }
    }

    if (fillFlag_local.isSet(Mesh::FILL_BOUND)) {
      for (int i = 0; i < 3; i++)
	boundary[10 + i] = elInfoOld->getBoundary(10 + cv[i]);
      
      boundary[13] = elInfoOld->getBoundary(4);
      
      boundary[0] = INTERIOR;
      boundary[1] = elInfoOld->getBoundary(cv[1]);
      boundary[2] = elInfoOld->getBoundary(cv[2]);
      boundary[3] = elInfoOld->getBoundary(ochild);
      
      int geoFace = mesh->getGeo(FACE);

      ce = const_cast<int*>(Tetrahedron::childEdge[el_type_local][ichild]);
      for (int iedge = 0; iedge < 4; iedge++)
	boundary[geoFace + iedge] = elInfoOld->getBoundary(geoFace + ce[iedge]);      
      for (int iedge = 4; iedge < 6; iedge++)
	boundary[geoFace + iedge] = elInfoOld->getBoundary(5 - cv[iedge - 3]);

      if (elInfoOld->getProjection(0) &&
	  elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
	
	projection[0] = elInfoOld->getProjection(0);      
      } else { // boundary projection
	projection[0] = NULL;
	projection[1] = elInfoOld->getProjection(cv[1]);
	projection[2] = elInfoOld->getProjection(cv[2]);
	projection[3] = elInfoOld->getProjection(ochild);
	
	for (int iedge = 0; iedge < 4; iedge++)
	  projection[geoFace + iedge] = elInfoOld->getProjection(geoFace + ce[iedge]);
	for (int iedge = 4; iedge < 6; iedge++)
	  projection[geoFace + iedge] = elInfoOld->getProjection(5 - cv[iedge - 3]);
      }
    }

    
    if (fillFlag.isSet(Mesh::FILL_ORIENTATION)) {
      orientation = 
	(dynamic_cast<ElInfo3d*>(const_cast<ElInfo*>(elInfoOld)))->orientation 
	* Tetrahedron::childOrientation[el_type_local][ichild];
    }
  }


  mtl::dense2D<double>& ElInfo3d::getSubElemCoordsMat(int degree) const
  {
    FUNCNAME("ElInfo3d::getSubElemCoordsMat()");

    ERROR_EXIT("Not yet implemented!\n");

    return subElemMatrices[degree][refinementPath];
  }


  mtl::dense2D<double>& ElInfo3d::getSubElemGradCoordsMat(int degree) const
  {
    FUNCNAME("ElInfo3d::getSubElemGradCoordsMat()");

    ERROR_EXIT("Not yet implemented!\n");

    return subElemGradMatrices[degree][refinementPath];
  }

}