-
Thomas Witkowski authoredThomas Witkowski authored
ElInfo3d.cc 17.81 KiB
#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;
el_type = 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();
el_type = (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++) {
int i = 5 - cv[iedge - 3]; /* old vertex opposite new edge */
boundary_[geoFace + iedge] = elInfoOld->getBoundary(i);
}
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];
}
}
}