#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<MacroElement*>( mel);
element_ = const_cast<Element*>( 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_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
WorldVector<double> oppC;
int neighbours = mesh_->getGeo(NEIGH);
for (int i = 0; i < neighbours; i++) {
nb = NULL;
if ((mnb = const_cast<MacroElement*>( mel->getNeighbour(i)))) {
if (fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
oppC = mnb->coord[i];
}
nb = const_cast<Element*>( 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<Element*>( 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<WorldVector<double> >& grd_lam)
{
FUNCNAME("ElInfo1d::calcGrdLambda()");
testFlag(Mesh::FILL_COORDS);
WorldVector<double> 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<double>& x,
DimVec<double>* lambda) const
{
FUNCNAME("ElInfo1d::worldToCoord()");
double lmin;
double a = coord_[0][0];
double length = (coord_[1][0] - a);
int dim = mesh_->getDim();
static DimVec<double> 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<double> &normal)
{
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<double> &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 *elInfoOld)
{
FUNCNAME("ElInfo1d::fillElInfo()");
Element *nb;
Element *elem = elInfoOld->element_;
TEST_EXIT_DBG(elem->getChild(0))("no children?\n");
element_ = const_cast<Element*>(elem->getChild(ichild));
TEST_EXIT_DBG(element_)("missing child %d?\n", ichild);
macroElement_ = elInfoOld->macroElement_;
fillFlag_ = elInfoOld->fillFlag_;
parent_ = elem;
level = elInfoOld->level + 1;
iChild = ichild;
int neighbours = mesh_->getGeo(NEIGH);
if (fillFlag_.isSet(Mesh::FILL_COORDS) || fillFlag_.isSet(Mesh::FILL_DET) ||
fillFlag_.isSet(Mesh::FILL_GRD_LAMBDA)) {
const FixVec<WorldVector<double>, VERTEX> *old_coord = &(elInfoOld->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_NEIGH) || fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
WorldVector<double> 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<Element*>( elem->getChild(1-ichild));
if ( fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
oppC = elInfoOld->coord_[i];
}
} else {
nb = const_cast<Element*>( elInfoOld->getNeighbour(i));
if (nb && fillFlag_.isSet(Mesh::FILL_OPP_COORDS)) {
oppC = elInfoOld->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<Element*>( 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] = elInfoOld->getBoundary(ichild);
boundary_[1 - ichild] = INTERIOR;
if (elInfoOld->getProjection(0) &&
elInfoOld->getProjection(0)->getType() == VOLUME_PROJECTION) {
projection_[0] = elInfoOld->getProjection(0);
}
}
return;
}
void ElInfo1d::getRefSimplexCoords(VectorOfFixVecs<DimVec<double> > *coords) const
{
(*coords)[0][0] = 1.0;
(*coords)[0][1] = 0.0;
(*coords)[1][0] = 0.0;
(*coords)[1][1] = 1.0;
}
void ElInfo1d::getSubElementCoords(VectorOfFixVecs<DimVec<double> > *coords,
int iChild) const
{
if (iChild == 0) {
(*coords)[1][0] = ((*coords)[0][0] + (*coords)[1][0]) * 0.5;
(*coords)[1][1] = ((*coords)[0][1] + (*coords)[1][1]) * 0.5;
} else {
(*coords)[0][0] = ((*coords)[0][0] + (*coords)[1][0]) * 0.5;
(*coords)[0][1] = ((*coords)[0][1] + (*coords)[1][1]) * 0.5;
}
}
}