#include "ElInfo1d.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 ElInfo1d::fillMacroInfo(const MacroElement * mel) { FUNCNAME("ElInfo1d::fillMacroInfo()"); Element *nb; MacroElement *mnb; macroElement_ = const_cast( mel); element_ = const_cast( mel->getElement()); parent_ = NULL; level_ = 0; 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]; } } // if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) { // det = elGrdLambda(*Lambda); // } if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { WorldVector oppC; int neighbours = mesh_->getGeo(NEIGH); for (int i = 0; i < neighbours; i++) { nb = NULL; if ((mnb = const_cast( mel->getNeighbour(i)))) { if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { oppC = mnb->coord[i]; } nb = const_cast( mnb->getElement()); while (!(nb->isLeaf())) { // make nb nearest element if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { if (nb->getNewCoord(-1)) { oppC = *(nb->getNewCoord()); } else { oppC = (mel->coord[i] + oppC) * 0.5; } } nb = const_cast( nb->getChild(1-i)); } if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { oppCoord_[i] = oppC; } } neighbour_[i] = nb; oppVertex_[i] = nb ? i : -1; } } if (fillFlag_.isSet(Mesh::FILL_BOUND) ) { for (int i = 0; i < vertices; i++) boundary_[i] = mel->getBoundary(i); for (int i = 0; i < element_->getGeo(PROJECTION); i++) { projection_[i] = mel->getProjection(i); } } } /****************************************************************************/ /* compute gradients of basis functions on element; return the absulute */ /* value of the determinante from the transformation to the reference */ /* element */ /****************************************************************************/ double ElInfo1d::calcGrdLambda(DimVec >& grd_lam) const { FUNCNAME("ElInfo1d::calcGrdLambda()"); testFlag(Mesh::FILL_COORDS); WorldVector e; e = coord_[1]; e -= coord_[0]; double adet2 = e * e; if (adet2 < 1.0E-15) { MSG("det*det = %lf\n", adet2); for (int i = 0; i <= 1; i++) grd_lam[i] = 0.0; } else { grd_lam[1] = e * (1.0 / adet2); grd_lam[0] = grd_lam[1] * (-1.0); } return sqrt(adet2); } const int ElInfo1d::worldToCoord(const WorldVector& x, DimVec* lambda) const { FUNCNAME("ElInfo1d::worldToCoord()"); double lmin; double a = coord_[0][0]; double length = (coord_[1][0] - a); int dim = mesh_->getDim(); static DimVec vec(dim, NO_INIT); TEST_EXIT_DBG(lambda)("lambda must not be NULL\n"); TEST_EXIT_DBG(dim == 1)("dim!=1\n"); TEST_EXIT_DBG(dimOfWorld == dim)("not yet for DIM != DIM_OF_WORLD\n"); if (abs(length) < DBL_TOL) { ERROR_EXIT("length = %le; abort\n", length); return 0; } (*lambda)[1] = (x[0] - a) / length; (*lambda)[0] = 1.0 - (*lambda)[1]; int k = -1; 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; } /****************************************************************************/ /* calculate a facenormal of edge side of a triangle with coordinates */ /* coord; return the absulute value of the determinant from the */ /* transformation to the reference element */ /****************************************************************************/ double ElInfo1d::getNormal(int side, WorldVector &normal) const { normal = coord_[side] - coord_[(side + 1) % 2]; double det = norm(&normal); TEST_EXIT_DBG(det > 1.e-30)("det = 0 on side %d\n", side); normal *= 1.0 / det; return(det); } /****************************************************************************/ /* calculate the normal of the element for dim of world = 2 */ /* return the absulute value of the determinant from the */ /* transformation to the reference element */ /****************************************************************************/ double ElInfo1d::getElementNormal( WorldVector &elementNormal) const { FUNCNAME("ElInfo::getElementNormal()"); TEST_EXIT_DBG(dimOfWorld == 2) (" element normal only well defined for DIM_OF_WORLD = DIM + 1 !!"); elementNormal[0] = coord_[1][1] - coord_[0][1]; elementNormal[1] = coord_[0][0] - coord_[1][0]; double det = norm(&elementNormal); TEST_EXIT_DBG(det > 1.e-30)("det = 0"); elementNormal *= 1.0 / det; return(det); } void ElInfo1d::fillElInfo(int ichild, const ElInfo *elinfo_old) { FUNCNAME("ElInfo1d::fillElInfo()"); Element *nb; Element *elem = elinfo_old->element_; TEST_EXIT_DBG(elem->getChild(0))("no children?\n"); element_ = const_cast(elem->getChild(ichild)); TEST_EXIT_DBG(element_)("missing child %d?\n", ichild); macroElement_ = elinfo_old->macroElement_; fillFlag_ = elinfo_old->fillFlag_; parent_ = elem; level_ = elinfo_old->level_ + 1; int neighbours = mesh_->getGeo(NEIGH); if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) { const FixVec, VERTEX> *old_coord = &(elinfo_old->coord_); coord_[ichild] = (*old_coord)[ichild]; if (elem->getNewCoord(-1)) { coord_[1 - ichild] = *(elem->getNewCoord()); } else { coord_[1 - ichild] = ((*old_coord)[0] + (*old_coord)[1]) * 0.5; } } // if(fillFlag_.isSet(Mesh::FILL_DET) || fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) { // det = calcGrdLambda(*Lambda); // } if (fillFlag_.isSet(Mesh::FILL_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { WorldVector oppC; TEST_EXIT_DBG(fillFlag_.isSet(Mesh::FILL_COORDS)) ("FILL_OPP_COORDS only with FILL_COORDS\n"); for (int i = 0; i < neighbours; i++) { if (i != ichild) { nb = const_cast( elem->getChild(1-ichild)); if ( fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { oppC = elinfo_old->coord_[i]; } } else { nb = const_cast( elinfo_old->getNeighbour(i)); if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { oppC = elinfo_old->oppCoord_[i]; } } if (nb) { while (nb->getChild(0)) { // make nb nearest element if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { if (nb->getNewCoord(-1)) { oppC = *(nb->getNewCoord()); } else { oppC = (coord_[i] + oppC) * 0.5; } } nb = const_cast( nb->getChild(1-i)); } if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) { oppCoord_[i] = oppC; } } neighbour_[i] = nb; oppVertex_[i] = nb ? i : -1; } } if (fillFlag_.isSet(Mesh::FILL_BOUND)) { boundary_[ichild] = elinfo_old->getBoundary(ichild); boundary_[1-ichild] = INTERIOR; if (elinfo_old->getProjection(0) && elinfo_old->getProjection(0)->getType() == VOLUME_PROJECTION) { projection_[0] = elinfo_old->getProjection(0); } } return; } }