/*
XLiFE++ is an extended library of finite elements written in C++
Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
/*!
\file GeomRefPrism.cpp
\authors D. Martin, N.Kielbasiewicz
\since 15 dec 2002
\date 19 sep 2012
\brief Implementation of xlifepp::GeomRefPrism class members and related functions
*/
#include "GeomRefPrism.hpp"
#include "utils.h"
namespace xlifepp
{
//--------------------------------------------------------------------------------
// GeomRefPrism constructor for Geometric Reference Element
//--------------------------------------------------------------------------------
GeomRefPrism::GeomRefPrism()
: GeomRefElement(_prism, 0.5, over3_, 6, 9)
//! using 3d base constructor with shape, volume, centroid coords, number of vertices, edges
{
trace_p->push("GeomRefPrism::GeomRefPrism");
centroid_[2] = 0.5;
// coordinates of vertices
std::vector::iterator it_v(vertices_.begin());
vertex(it_v, 1., 0., 0.);
vertex(it_v, 0., 1., 0.);
vertex(it_v, 0., 0., 0.);
vertex(it_v, 1., 0., 1.);
vertex(it_v, 0., 1., 1.);
vertex(it_v, 0., 0., 1.);
// vertex numbering on edges
sideOfSideNumbering();
// vertex numbering and oriented edge numbering on faces
sideNumbering();
trace_p->pop();
}
/*
--------------------------------------------------------------------------------
sideOfSideNumbering
defines Geometric Reference Element local numbering of edge vertices
sideNumbering
defines Geometric Reference Element local numbering of vertices on faces
and of edges on faces
--------------------------------------------------------------------------------
*/
void GeomRefPrism::sideOfSideNumbering()
{
sideOfSideVertexNumbers_[0].push_back(1);
sideOfSideVertexNumbers_[0].push_back(2);
sideOfSideVertexNumbers_[1].push_back(2);
sideOfSideVertexNumbers_[1].push_back(3);
sideOfSideVertexNumbers_[2].push_back(3);
sideOfSideVertexNumbers_[2].push_back(1);
sideOfSideVertexNumbers_[3].push_back(4);
sideOfSideVertexNumbers_[3].push_back(5);
sideOfSideVertexNumbers_[4].push_back(5);
sideOfSideVertexNumbers_[4].push_back(6);
sideOfSideVertexNumbers_[5].push_back(6);
sideOfSideVertexNumbers_[5].push_back(4);
sideOfSideVertexNumbers_[6].push_back(1);
sideOfSideVertexNumbers_[6].push_back(4);
sideOfSideVertexNumbers_[7].push_back(2);
sideOfSideVertexNumbers_[7].push_back(5);
sideOfSideVertexNumbers_[8].push_back(3);
sideOfSideVertexNumbers_[8].push_back(6);
sideOfSideNumbers_[0].push_back(1);
sideOfSideNumbers_[0].push_back(-3);
sideOfSideNumbers_[0].push_back(-2);
sideOfSideNumbers_[1].push_back(1);
sideOfSideNumbers_[1].push_back(7);
sideOfSideNumbers_[1].push_back(4);
sideOfSideNumbers_[1].push_back(-8);
sideOfSideNumbers_[2].push_back(-9);
sideOfSideNumbers_[2].push_back(-5);
sideOfSideNumbers_[2].push_back(-8);
sideOfSideNumbers_[2].push_back(2);
sideOfSideNumbers_[3].push_back(-3);
sideOfSideNumbers_[3].push_back(-7);
sideOfSideNumbers_[3].push_back(6);
sideOfSideNumbers_[3].push_back(9);
sideOfSideNumbers_[4].push_back(4);
sideOfSideNumbers_[4].push_back(-5);
sideOfSideNumbers_[4].push_back(-6);
}
void GeomRefPrism::sideNumbering()
{
for (number_t i = 1; i < nbSides_ - 1; i++)
{
sideShapeTypes_[i] = _quadrangle;
}
sideShapeTypes_[0] = _triangle;
sideShapeTypes_[nbSides_ - 1] = _triangle;
sideVertexNumbers_[0].push_back(1);
sideVertexNumbers_[0].push_back(3);
sideVertexNumbers_[0].push_back(2);
sideVertexNumbers_[1].push_back(1);
sideVertexNumbers_[1].push_back(2);
sideVertexNumbers_[1].push_back(5);
sideVertexNumbers_[1].push_back(4);
sideVertexNumbers_[2].push_back(2);
sideVertexNumbers_[2].push_back(3);
sideVertexNumbers_[2].push_back(6);
sideVertexNumbers_[2].push_back(5);
sideVertexNumbers_[3].push_back(3);
sideVertexNumbers_[3].push_back(1);
sideVertexNumbers_[3].push_back(4);
sideVertexNumbers_[3].push_back(6);
sideVertexNumbers_[4].push_back(4);
sideVertexNumbers_[4].push_back(5);
sideVertexNumbers_[4].push_back(6);
}
//! Returns edge length, face area or element volume
real_t GeomRefPrism::measure(const dimen_t d, const number_t sideNum) const
{
real_t ms = 0.;
switch (d)
{
case 0: ms = 1.; break;
case 1:
switch (sideNum)
{
case 1: case 4: ms = sqrtOf2_; break;
case 2: case 3: case 5: case 6: case 7: case 8: case 9: ms = 1.; break;
default: noSuchSideOfSide(sideNum);
}
break;
case 2:
switch (sideNum)
{
case 1: case 5: ms = .5; break;
case 3: case 4: ms = 1.; break;
case 2: ms = sqrtOf2_; break;
default: noSuchSide(sideNum);
}
break;
case 3: ms = measure_; break;
default: break;
}
return ms;
}
// /*
// --------------------------------------------------------------------------------
// Returns a quadrature rule built on an edge from a 1d quadrature formula
// or a quadrature rule built on an face from a triangle or quadrangle quadrature formula
// --------------------------------------------------------------------------------
// */
// void GeomRefPrism::sideQuadrature(const QuadratureRule& q1, QuadratureRule& qr, const number_t sideNum, const dimen_t d) const
// {
// /*
// q1 : input 1d or 2d (triangle or quadrangle) quadrature rule
// q2 : ouput 3d quadrature rule set on the edge or face
// sideNum : local number of edge (sideNum = 1,2,...,9) or face (sideNum = 1,2,3,4,5)
// d : side dimension (d = 1 for an edge, d = 2 for a face)
// */
// std::vector::const_iterator c_1(q1.point()), w_1(q1.weight());
// std::vector::iterator c_i(qr.point()), w_i(qr.weight());
// switch ( d )
// {
// case 2:
// switch ( sideNum )
// {
// case 1: // face (x3 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
// qr.point(c_i, *c_1, *(c_1+1), 0., w_i, *w_1);
// break;
// case 2: // face (x1 + x2 = 1)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
// qr.point(c_i, *c_1, 1.-*c_1, *(c_1+1), w_i, *w_1*sqrt_of_2_);
// break;
// case 3: // face (x1 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
// qr.point(c_i, 0., *c_1, *(c_1+1), w_i, *w_1);
// break;
// case 4: // face (x2 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
// qr.point(c_i, *c_1, 0., *(c_1+1), w_i, *w_1);
// break;
// case 5: // face (x3 = 1)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
// qr.point(c_i, *c_1, *(c_1+1), 1., w_i, *w_1);
// break;
// default: noSuchFace(sideNum); break;
// }
// break;
// case 1:
// switch ( sideNum )
// {
// // edges in plane x3 = 0
// case 1: // edge (x1 + x2 = 1 & x3 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, *c_1, 1.-*c_1, 0., w_i, *w_1*sqrt_of_2_);
// break;
// case 2: // edge (x1 = 0 & x3 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, 0., *c_1, 0., w_i, *w_1);
// break;
// case 3: // edge (x2 = 0 & x3 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, *c_1, 0., 0., w_i, *w_1);
// break;
// // edges in plane x3 = 1
// case 4: // edge (x1 + x2 = 1 & x3 = 1)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, *c_1, 1.-*c_1, 1., w_i, *w_1*sqrt_of_2_);
// break;
// case 5: // edge (x1 = 0 & x3 = 1)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, 0., *c_1, 1., w_i, *w_1);
// break;
// case 6: // edge (x2 = 0 & x3 = 1)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, *c_1, 0., 1., w_i, *w_1);
// break;
// // edges orthogonal to x3 = 0 plane
// case 7: // edge (x1 = 1 & x2 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, 1., 0., *c_1, w_i, *w_1);
// break;
// case 8: // edge (x1 = 0 & x2 = 1)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, 0., 1., *c_1, w_i, *w_1);
// break;
// case 9: // edge (x1 = 0 & x2 = 0)
// for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
// qr.point(c_i, 0., 0., *c_1, w_i, *w_1);
// break;
// default: noSuchEdge(sideNum);
// break;
// }
// default: break;
// }
// }
//! returns tangent vector(s) on edge sideNum (=1,2,...,9) or face sideNum (=1,2,3,4,5)
void GeomRefPrism::tangentVector(const std::vector& jacobianMatrix, std::vector >& tgv, const number_t sideNum, const dimen_t d) const
{
std::vector::iterator it_tgv0(tgv[0].begin()), it_tgv1(it_tgv0);
if ( d > 1) { it_tgv1 = tgv[1].begin(); }
// std::vector::const_iterator it_jm(jacobianMatrix.begin());
noSuchFunction("tangentVector");
// switch ( d )
// {
// case 2:
// switch ( sideNum )
// {
// case 1: //
// break;
// case 2: //
// break;
// case 3: //
// break;
// case 4: //
// break;
// case 5: //
// break;
// default: noSuchSide(sideNum); break;
// }
// case 1:
// switch ( sideNum )
// {
// case 1: //
// break;
// case 2: //
// break;
// case 3: //
// break;
// case 4: //
// break;
// case 5: //
// break;
// case 6: //
// break;
// case 7: //
// break;
// case 8: //
// break;
// case 9: //
// break;
// default: noSuchSideOfSide(sideNum); break;
// }
// break;
// default:
// break;
// }
}
//! returns local number of edge bearing vertices with local numbers v1 and v2
number_t GeomRefPrism::sideWithVertices(const number_t vn1, const number_t vn2) const
{
if(vn1==vn2) noSuchSide(vn1,vn2);
number_t v1=vn1, v2=vn2;
if(v1>v2) {v1=vn2; v2=vn1;}
switch(v1)
{
case 1 :
{
if(v2==2) return 1;
if(v2==3) return 3;
if(v2==4) return 7;
}
break;
case 2 :
{
if(v2==3) return 2;
if(v2==5) return 8;
}
break;
case 3 : if(v2==6) return 9; break;
case 4 :
{
if(v2==5) return 4;
if(v2==6) return 6;
}
break;
case 5 : if(v2==6) return 5;
}
noSuchSide(vn1,vn2);
return 0;
}
/*! returns local number of face bearing vertices with local numbers v1, v2, v3 and v4 if not null
1: 1 2 3 100v1+10v2+v3 123
2: 1 2 5 4 124 125 145 245
3: 2 3 6 5 235 236 256 356
4: 1 3 6 4 134 136 146 346
5: 4 5 6 456
*/
number_t GeomRefPrism::sideWithVertices(const number_t vn1, const number_t vn2, const number_t vn3, const number_t vn4) const
{
std::set vs; vs.insert(vn1); vs.insert(vn2); vs.insert(vn3);
if(vn4>0)
{
vs.insert(vn4);
if(vs.size()!=4 || *vs.begin()<1 || *vs.rbegin()>6) noSuchSide(vn1,vn2,vn3,vn4);
}
else if(vs.size()!=3 || *vs.begin()<1 || *vs.rbegin()>6) noSuchSide(vn1,vn2,vn3);
std::map hm;
std::map::iterator itm;
number_t s=0;
std::set::iterator itv=vs.begin();
if(vn4==0)
{
hm[123]=1; hm[124]=2; hm[125]=2; hm[145]=2; hm[245]=2; hm[235]=3; hm[236]=3; hm[256]=3; hm[356]=3;
hm[134]=4; hm[136]=4; hm[146]=4; hm[346]=4; hm[456]=5;
s=100**itv++; s+=10**itv++; s+=*itv;
itm=hm.find(s);
if(itm!=hm.end()) return itm->second; else noSuchSide(vn1,vn2,vn3);
}
hm[1245]=2; hm[2356]=3; hm[1346]=4;
s=1000**itv++; s+=100**itv++; s+=10**itv++;s+=*itv;
itm=hm.find(s);
if(itm!=hm.end()) return itm->second; else noSuchSide(vn1,vn2,vn3,vn4);
return 0;
}
//! test if a point belongs to current element
bool GeomRefPrism::contains(std::vector& p, real_t tol) const
{
real_t x=p[0], y=p[1], z=p[2];
return (x >= -tol) && (x <= 1.+tol) && (y >= -tol) && (y <= 1.+tol) && (z >= -tol) && (z <= 1.+tol)
&& (x+y<=1.+tol);
}
} // end of namespace xlifepp