19 static LoggerPtr logger(Logger::getLogger(
"pz.topology.pztriangle"));
28 T qsi = loc[0], eta = loc[1];
29 phi(0,0) = 1.0-qsi-eta;
33 dphi(0,0) = dphi(1,0) = -1.0 + qsi;
34 dphi(0,1) = dphi(1,2) = 1.0 + qsi;
35 dphi(1,1) = dphi(0,2) = qsi;
38 int TPZTriangle::SideNodes[3][2] = { {0,1},{1,2},{2,0} };
39 int TPZTriangle::FaceNodes[1][3] = { {0,1,2} };
43 REAL TPZTriangle::gTrans2dT[6][2][2] = {
44 { { 1., 0.},{ 0., 1.} },
45 { { 0., 1.},{ 1., 0.} },
46 { { 0., 1.},{-1.,-1.} },
47 { {-1.,-1.},{ 0., 1.} },
48 { {-1.,-1.},{ 1., 0.} },
49 { { 1., 0.},{-1.,-1.} }
53 void TPZTriangle::BlendFactorForSide(
const int &side,
const TPZVec<T> &xi, T &blendFactor,
56 blendFactorDxi.
Resize(TPZTriangle::Dimension, (T) 0);
58 std::ostringstream sout;
59 if(side < NCornerNodes || side >= NSides){
60 sout<<
"The side\t"<<side<<
"is invalid. Aborting..."<<std::endl;
64 sout<<
"The method BlendFactorForSide expects the point xi to correspond to coordinates of a point";
65 sout<<
" inside the parametric domain. Aborting...";
68 if(!CheckProjectionForSingularity(side,xi)){
69 sout<<
"The projection of xi "<<xi[0]<<
" "<<xi[1]<<
" to side "<<side<<
" is singular."<<std::endl;
70 sout<<
"This should have been caught by MapToSide method. Aborting..."<<std::endl;
72 if(!sout.str().empty()){
73 PZError<<std::endl<<sout.str()<<std::endl;
81 if(!CheckProjectionForSingularity(side,xi)){
82 std::cout<<
"Side projection is not regular and it should have been checked earlier. Aborting.."<<std::endl;
85 for(
int i = 0; i < blendFactorDxi.
size(); i++) blendFactorDxi[i] = 0;
91 TPZTriangle::TShape(xi,phi,dphi);
112 blendFactor = phi(i,0) + phi((i+1)%NCornerNodes,0);
113 blendFactor *= blendFactor;
114 blendFactorDxi[0] = 2 * ( phi(i,0) + phi((i+1)%NCornerNodes,0) ) * ( dphi(0,i) + dphi(0,(i+1)%NCornerNodes) );
115 blendFactorDxi[1] = 2 * ( phi(i,0) + phi((i+1)%NCornerNodes,0) ) * ( dphi(1,i) + dphi(1,(i+1)%NCornerNodes) );
135 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
136 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
137 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,-99}}
140 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
141 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
142 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,0,-99}}
145 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
146 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
147 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,1,-99}}
150 {{0.5,0,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0,-99}}
153 {{-0.5,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0.5,-99}}
156 {{0,-0.5,-99},{-99,-99,-99},{-99,-99,-99},{0,0.5,-99}}
159 {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-99,-99,-99}}
164 {.0,0.},{1.0,.0},{0.,1.0},
165 {.5,0.},{0.5,.5},{0.,0.5},
207 static int vectorsideorder [14] = {0,1,3,1,2,4,2,0,5,3,4,5,6,6};
210 static int bilinearounao [14] = {0,0,0,0,0,0,0,0,0,1,1,1,1,1};
212 static int direcaoksioueta [14] = {0,0,0,0,0,0,0,0,0,0,0,0,0,1};
214 int TPZTriangle::fPermutations [6][7] =
224 REAL TPZTriangle::fTangentVectors [12][2] =
240 int TPZTriangle::NBilinearSides()
247 int TPZTriangle::NSideNodes(
int side)
249 return nsidenodes[side];
253 const REAL qsi = pt[0];
254 const REAL eta = pt[1];
255 if( ( qsi <= 1. + tol ) && ( qsi >= 0. - tol ) &&
256 ( eta <= 1. + tol ) && ( eta >= 0. - tol ) &&
257 ( eta <= 1. - qsi + tol ) ){
266 bool TPZTriangle::CheckProjectionForSingularity(
const int &side,
const TPZVec<T> &xiInterior) {
269 T qsi = xiInterior[0]; T eta = xiInterior[1];
278 if(
fabs((T)(eta - 1.)) < zero)
return false;
280 if((T)(qsi+eta) < (T)zero)
return false;
282 if(
fabs((T)(qsi - 1.)) < zero)
return false;
287 cout <<
"Cant compute CheckProjectionForSingularity method in TPZTriangle class!\nParameter (SIDE) must be 3, 4 or 5!\nMethod Aborted!\n";
296 T qsi = InternalPar[0]; T eta = InternalPar[1];
299 if(!CheckProjectionForSingularity(side,InternalPar)){
300 std::cout<<
"Side projection is not regular and it should have been checked earlier. Aborting.."<<std::endl;
312 SidePar[0] = 2.*qsi/(1.-eta) - 1.;
313 JacToSide(0,0) = 2./(1.-eta); JacToSide(0,1) = 2.*qsi/((1.-eta)*(1.-eta));
317 SidePar[0] = 1. - 2.*qsi/(qsi + eta);
318 JacToSide(0,0) = -2.*eta/((qsi+eta)*(qsi+eta)); JacToSide(0,1) = 2.*qsi/((qsi+eta)*(qsi+eta));
322 SidePar[0] = 1. - 2.*eta/(1.-qsi);
323 JacToSide(0,0) = -2.*eta/((1.-qsi)*(1.-qsi)); JacToSide(0,1) = -2./(1.-qsi);
326 SidePar = InternalPar;
332 void TPZTriangle::ParametricDomainNodeCoord(
int node,
TPZVec<REAL> &nodeCoord)
334 if(node > NCornerNodes)
338 nodeCoord.
Resize(Dimension, 0.);
369 REAL
val = (REAL) rand() / (RAND_MAX);
371 val = (1.-pt[0]) * (REAL) rand() / (RAND_MAX);
378 if(sidefrom <0 || sidefrom >= NSides || sideto <0 || sideto >= NSides) {
379 PZError <<
"TPZTriangle::HigherDimensionSides sidefrom "<< sidefrom <<
380 ' ' << sideto << endl;
383 if(sidefrom == sideto) {
386 if(sidefrom == NSides-1) {
387 return TransformElementToSide(sideto);
389 if (sideto == NSides-1) {
390 return TransformSideToElement(sidefrom);
393 int nhigh = nhighdimsides[sidefrom];
395 for(is=0; is<nhigh; is++) {
396 if(highsides[sidefrom][is] == sideto) {
397 int dfr = sidedimension[sidefrom];
398 int dto = sidedimension[sideto];
401 for(i=0; i<dto; i++) {
402 for(j=0; j<dfr; j++) {
403 trans.
Mult()(i,j) = sidetosidetransforms[sidefrom][is][j][i];
405 trans.
Sum()(i,0) = sidetosidetransforms[sidefrom][is][3][i];
410 PZError <<
"TPZTriangle::SideToSideTransform highside not found sidefrom " 411 << sidefrom <<
' ' << sideto << endl;
415 int TPZTriangle::SideNodeLocId(
int side,
int node)
417 if(side<3 && node == 0)
return side;
418 if(side>=3 && side<6 && node <2)
return (side-3+node) %3;
419 if(side==6 && node <3)
return node;
420 PZError <<
"TPZTriangle::SideNodeLocId inconsistent side or node " << side
421 <<
' ' << node << endl;
428 int nsidecon = NContainedSides(side);
430 for(is=0; is<nsidecon-1; is++)
431 smallsides.
Push(ContainedSideLocId(side,is));
434 void TPZTriangle::LowerDimensionSides(
int side,
TPZStack<int> &smallsides,
int DimTarget)
437 int nsidecon = NContainedSides(side);
438 for(
int is = 0; is < nsidecon - 1; is++) {
439 if (SideDimension(ContainedSideLocId(side,is)) == DimTarget) smallsides.
Push(ContainedSideLocId(side,is));
445 if(side <0 || side >= NSides) {
446 PZError <<
"TPZTriangle::HigherDimensionSides side "<< side << endl;
449 for(is=0; is<nhighdimsides[side]; is++) high.
Push(highsides[side][is]);
455 if(dimension<0 || dimension> 2) {
456 PZError <<
"TPZTriangle::NumSides. Bad parameter i.\n";
459 if(dimension==0)
return 3;
460 if(dimension==1)
return 3;
461 if(dimension==2)
return 1;
464 int TPZTriangle::SideDimension(
int side) {
465 if(side<0 || side >= NSides) {
466 PZError <<
"TPZTriangle::SideDimension side " << side << endl;
469 return sidedimension[side];
473 if (center.
size()!=Dimension) {
477 for(i=0; i<Dimension; i++) {
478 center[i] = MidSideNode[side][i];
484 if(side<0 || side>6){
485 PZError <<
"TPZTriangle::TransformElementToSide called with side error\n";
504 t.
Mult()(0,0) = -1.0;
508 t.
Mult()(0,0) = -1.0;
509 t.
Mult()(0,1) = -2.0;
525 if(side<0 || side>6){
526 PZError <<
"TPZTriangle::TransformSideToElement side out range\n";
547 t.
Mult()(0,0) = -0.5;
553 t.
Mult()(1,0) = -0.5;
564 TPZIntPoints * TPZTriangle::CreateSideIntegrationRule(
int side,
int order){
565 if(side < 0 || side>6) {
566 PZError <<
"TPZTriangle::CreateSideIntegrationRule wrong side = " << side <<
".\n";
570 if(side<6)
return new TPZInt1d(order);
600 int TPZTriangle::NumSides() {
606 int TPZTriangle::NContainedSides(
int side) {
607 if(side<0 || side>6) {
608 PZError <<
"TPZShapeTriang::NContainedSides. Bad parameter side = " << side <<
".\n";
618 int TPZTriangle::ContainedSideLocId(
int side,
int c) {
627 if(!c)
return side-3;
628 if(c==1)
return (side-2)%3;
629 if(c==2)
return side;
633 PZError <<
"TPZShapeTriang::ContainedSideLocId, connect = " << c <<
".\n";
647 id0 = (
id[0] <
id[1]) ? 0 : 1;
648 minid = (
id[2] <
id[id0]) ? 2 : id0;
649 id0 = (minid + 1) % 3;
650 id1 = (minid + 2) % 3;
652 if (
id[id0] <
id[id1]) {
654 if (minid == 0)
return 0;
655 if (minid == 1)
return 2;
656 if (minid == 2)
return 4;
661 if (minid == 0)
return 1;
662 if (minid == 1)
return 3;
663 if (minid == 2)
return 5;
686 int in1 = ContainedSideLocId(side,0);
687 int in2 = ContainedSideLocId(side,1);
688 return id[in1]<
id[in2] ? 0 : 1;
693 return GetTransformId(
id);
706 void TPZTriangle::GetHDivGatherPermute(
int transformid,
TPZVec<int> &permute)
709 if (permute.
size() != 7 || transformid >= 6 || transformid < 0) {
713 for (
int i=0; i<7; i++) {
714 permute[i] = fPermutations[transformid][i];
718 if (transformid%2 ==1) dir = -1;
719 int runsmall = transformid/2;
720 int runlarge = runsmall+3;
728 for (
int is=0; is<3; is++) {
729 permute[is] = runsmall;
730 permute[is+3] = runlarge;
733 if (dir == 1 && runsmall > 2) {
737 else if(dir == -1 && runsmall < 0)
741 if (dir == -1 && runlarge < 3) {
755 void TPZTriangle::GetSideHDivPermutation(
int transformationid,
TPZVec<int> &permgather)
759 for (
int i=0; i<7; i++)
761 permgather[i] = fPermutations[transformationid][i];
1060 template <
class TVar>
1066 for (
int i=0; i<3; i++) {
1069 vdiag[i] = (gradx(i,0)-gradx(i,1));
1086 NormalScales[0] = 2./Nv1;
1087 NormalScales[1] = 2./Nvdiag;
1088 NormalScales[2] = 2./Nv2;
1091 for (
int i=0; i<3; i++) {
1097 for (
int i=0; i<3; i++)
1099 directions(i,0) = -v2[i]*Nv1*NormalScales[0];
1100 directions(i,1) = (v1[i]-v2[i])*Nv1*NormalScales[0];
1101 directions(i,2) = (directions(i,0)+directions(i,1))/2.;
1102 directions(i,3) = v1[i]*Nvdiag*NormalScales[1];
1103 directions(i,4) = v2[i]*Nvdiag*NormalScales[1];
1104 directions(i,5) = (directions(i,3)+directions(i,4))/2.;
1105 directions(i,6) = (v2[i]-v1[i])*Nv2*NormalScales[2];
1106 directions(i,7) = -v1[i]*Nv2*NormalScales[2];
1107 directions(i,8) = (directions(i,6)+directions(i,7))/2.;
1109 directions(i,9) = v1[i]*Nv1*NormalScales[0];
1110 directions(i,10) = (v2[i]-v1[i])*Nvdiag*NormalScales[1];
1111 directions(i,11) = -v2[i]*Nv2*NormalScales[2];
1113 directions(i,12) = v1[i]*Nv2*NormalScales[0];
1114 directions(i,13) = v2[i]*Nv1*NormalScales[1];
1118 template <
class TVar>
1121 const auto dim = gradx.
Rows();
1123 for (
int i=0; i<dim; i++) {
1127 constexpr
auto nEdges{3};
1128 constexpr REAL faceArea{TPZTriangle::RefElVolume()};
1129 constexpr REAL edgeLength[nEdges]{1,M_SQRT2,1};
1130 const int facePermute = transformationIds[nEdges];
1132 for(
auto iEdge = 0; iEdge < nEdges; iEdge++){
1133 edgeSign[iEdge] = transformationIds[iEdge] == 0 ? 1 : -1;
1135 for (
int i=0; i<dim; i++)
1140 directions(i,0) = (v1[i]) * edgeSign[0] / edgeLength[0];
1141 directions(i,1) = (v1[i]+v2[i]) * edgeSign[0] / edgeLength[0];
1142 directions(i,2) = (v2[i]*M_SQRT2) * edgeSign[1] / edgeLength[1];
1143 directions(i,3) = (-v1[i]*M_SQRT2) * edgeSign[1] / edgeLength[1];
1144 directions(i,4) = (v1[i] + v2[i]) * -1 * edgeSign[2] / edgeLength[2];
1145 directions(i,5) = (-v2[i]) * edgeSign[2] / edgeLength[2];
1148 directions(i,6) = (v1[i]) * edgeSign[0] / edgeLength[0];
1149 directions(i,7) = ( (v2[i] - v1[i]) * M_SQRT1_2) * edgeSign[1] / edgeLength[1];
1150 directions(i,8) = (-v2[i]) * edgeSign[2] / edgeLength[2];
1155 directions(i,9) = v2[i] * edgeSign[0] /faceArea;
1156 directions(i,10) = ((-v2[i]-v1[i]) * M_SQRT1_2) * edgeSign[1] /faceArea;
1157 directions(i,11) = v1[i] * edgeSign[2] /faceArea;
1160 directions(i,12) = fTangentVectors[2*facePermute][i];
1161 directions(i,13) = fTangentVectors[2*facePermute + 1][i];
1165 template <
class TVar>
1168 for (
auto i=0; i< Dimension; i++){
1170 v1[i] = fTangentVectors[2*transformationId][i];
1171 v2[i] = fTangentVectors[2*transformationId + 1][i];
1177 int nsides = NumSides()*2;
1183 for (
int is = 0; is<nsides; is++)
1185 sides[is] = vectorsideorder[is];
1186 dir[is] = direcaoksioueta[is];
1187 bilounao[is] = bilinearounao[is];
1193 int nsides = NumSides()*2;
1199 for (
int is = 0; is<nsides; is++)
1201 sides[is] = vectorsideorder[is];
1202 dir[is] = direcaoksioueta[is];
1203 bilounao[is] = bilinearounao[is];
1206 for (
int i=0; i<Dimension*NumSides(); i++) {
1207 sidevectors[i] = vectorsideorder[i];
1211 int TPZTriangle::ClassId()
const{
1212 return Hash(
"TPZTriangle");
1219 void TPZTriangle::Write(
TPZStream& buf,
int withclassid)
const {
1230 template bool pztopology::TPZTriangle::CheckProjectionForSingularity<REAL>(
const int &side,
const TPZVec<REAL> &xiInterior);
1234 template void pztopology::TPZTriangle::BlendFactorForSide<REAL>(
const int &,
const TPZVec<REAL> &, REAL &, TPZVec<REAL> &);
1242 template void pztopology::TPZTriangle::ComputeHCurlFaceDirections<REAL>(TPZVec<REAL> &v1, TPZVec<REAL> &v2,
int transformationId);
1245 template bool pztopology::TPZTriangle::CheckProjectionForSingularity<Fad<REAL> >(
const int &side,
const TPZVec<Fad<REAL> > &xiInterior);
1249 template void pztopology::TPZTriangle::BlendFactorForSide<Fad<REAL>>(
const int &,
const TPZVec<Fad<REAL>> &,
Fad<REAL> &,
1250 TPZVec<Fad<REAL>> &);
1257 template void pztopology::TPZTriangle::ComputeHCurlFaceDirections<Fad<REAL>>(TPZVec<Fad<REAL>> &v1, TPZVec<Fad<REAL>> &v2,
int transformationId);
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
static int direcaoksioueta[14]
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
Handles the numerical integration for two-dimensional problems using triangular elements. Numerical Integration.
Implements a vector class which allows to use external storage provided by the user. Utility.
static bool IsInParametricDomain(const TPZVec< REAL > &pt, REAL tol=pztopology::gTolerance)
Verifies if the parametric point pt is in the element parametric domain.
clarg::argInt dimension("-d", "Matrices dimension M x M", 1000)
static int bilinearounao[14]
virtual void Identity()
Converts the matrix in an identity matrix.
This class implements a simple vector storage scheme for a templated class T. Utility.
REAL val(STATE &number)
Returns value of the variable.
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Contains the TPZTriangle class which defines the topology of a triangle.
static REAL tTriang[14][2]
Handles the numerical integration for one-dimensional problems. Numerical Integration.
#define LOGPZ_FATAL(A, B)
Define log for fatal errors (cout)
Abstract class defining integration rules. Numerical Integration.
int Zero() override
Makes Zero all the elements.
int64_t size() const
Returns the number of elements of the vector.
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Groups all classes defining the structure of the master element.
void Push(const T object)
Pushes a copy of the object on the stack.
#define DebugStop()
Returns a message to user put a breakpoint in.
int64_t Rows() const
Returns number of rows.
static int sidedimension[7]
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
static REAL MidSideNode[7][3]
int32_t Hash(std::string str)
static int vectorsideorder[14]
static REAL bTriang[14][2]
static Tvar Norm(const TPZVec< Tvar > &vetor)
Returns the L2-norm of the vector.
Integration rule for one point. Numerical Integration.
MElementType
Define the element types.
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Defines the interface for saving and reading data. Persistency.
static int highsides[7][3]
static REAL sidetosidetransforms[7][3][4][3]
static int nhighdimsides[7]
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
#define PZError
Defines the output device to error messages and the DebugStop() function.