Select Git revision
Callbacks.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Christophe Geuzaine authoredChristophe Geuzaine authored
HierarchicalBasisH1Brick.cpp 22.37 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
//
// Contributed by Ismail Badia.
// Reference : "Higher-Order Finite Element Methods"; Pavel Solin, Karel
// Segeth ,
// Ivo Dolezel , Chapman and Hall/CRC; Edition : Har/Cdr (2003).
#include "HierarchicalBasisH1Brick.h"
HierarchicalBasisH1Brick::HierarchicalBasisH1Brick(int order)
{
_pb1 = order;
_pb2 = order;
_pb3 = order;
for(int i = 0; i < 12; i++) { _pOrderEdge[i] = order; }
for(int j = 0; j < 6; j++) {
_pOrderFace1[j] = order;
_pOrderFace2[j] = order;
}
_nvertex = 8;
_nedge = 12;
_nfaceQuad = 6;
_nfaceTri = 0;
_nVertexFunction = 8;
_nEdgeFunction = 12 * order - 12;
_nQuadFaceFunction = 6 * (order - 1) * (order - 1);
_nTriFaceFunction = 0;
_nBubbleFunction = (order - 1) * (order - 1) * (order - 1);
}
HierarchicalBasisH1Brick::~HierarchicalBasisH1Brick() {}
double HierarchicalBasisH1Brick::_affineCoordinate(const int &j,
const double &u,
const double &v,
const double &w)
{
switch(j) {
case(1): return 0.5 * (1 + u);
case(2): return 0.5 * (1 - u);
case(3): return 0.5 * (1 + v);
case(4): return 0.5 * (1 - v);
case(5): return 0.5 * (1 + w);
case(6): return 0.5 * (1 - w);
default: throw std::string("j must be : 1<=j<=6");
}
}
inline void HierarchicalBasisH1Brick::_someProduct(double const &u,
double const &v,
double const &w,
std::vector<double> &product,
std::vector<double> &lambda)
{
lambda[0] = _affineCoordinate(1, u, v, w);
lambda[1] = _affineCoordinate(2, u, v, w);
lambda[2] = _affineCoordinate(3, u, v, w);
lambda[3] = _affineCoordinate(4, u, v, w);
lambda[4] = _affineCoordinate(5, u, v, w);
lambda[5] = _affineCoordinate(6, u, v, w);
product[0] = lambda[3] * lambda[5];
product[1] = lambda[1] * lambda[5];
product[2] = lambda[1] * lambda[3];
product[3] = lambda[0] * lambda[5];
product[4] = lambda[3] * lambda[0];
product[5] = lambda[2] * lambda[5];
product[6] = lambda[2] * lambda[0];
product[7] = lambda[2] * lambda[1];
product[8] = lambda[3] * lambda[4];
product[9] = lambda[4] * lambda[1];
product[10] = lambda[4] * lambda[0];
product[11] = lambda[4] * lambda[2];
}
void HierarchicalBasisH1Brick::generateBasis(double const &u, double const &v,
double const &w,
std::vector<double> &vertexBasis,
std::vector<double> &edgeBasis,
std::vector<double> &faceBasis,
std::vector<double> &bubbleBasis)
{
std::vector<double> product(12, 0);
std::vector<double> lambda(6, 0);
HierarchicalBasisH1Brick::_someProduct(u, v, w, product, lambda);
// vertex shape functions:
vertexBasis[0] = lambda[1] * product[0];
vertexBasis[1] = lambda[0] * product[0];
vertexBasis[2] = lambda[0] * product[5];
vertexBasis[3] = lambda[1] * product[5];
vertexBasis[4] = lambda[1] * product[8];
vertexBasis[5] = lambda[0] * product[8];
vertexBasis[6] = lambda[0] * product[11];
vertexBasis[7] = lambda[1] * product[11];
std::vector<double> lkVectorU(_pb1 - 1);
std::vector<double> lkVectorV(_pb2 - 1);
std::vector<double> lkVectorW(_pb3 - 1);
for(int it = 2; it <= _pb1; it++) {
lkVectorU[it - 2] = OrthogonalPoly::EvalLobatto(it, u);
}
for(int it = 2; it <= _pb2; it++) {
lkVectorV[it - 2] = OrthogonalPoly::EvalLobatto(it, v);
}
for(int it = 2; it <= _pb3; it++) {
lkVectorW[it - 2] = OrthogonalPoly::EvalLobatto(it, w);
}
// edge shape functions:
int indexEdgeBasis = 0;
std::vector<double> *vectorTarget1(0);
for(int iEdge = 0; iEdge < _nedge; iEdge++) {
switch(iEdge) {
case(0):
case(5):
case(8):
case(11): vectorTarget1 = &lkVectorU; break;
case(1):
case(3):
case(9):
case(10): vectorTarget1 = &lkVectorV; break;
case(2):
case(4):
case(6):
case(7): vectorTarget1 = &lkVectorW; break;
}
for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] - 1;
indexEdgeFunc++) {
edgeBasis[indexEdgeBasis] =
(*vectorTarget1)[indexEdgeFunc] * product[iEdge];
indexEdgeBasis++;
}
}
// face shape functions:
int indexFaceFunction = 0;
std::vector<double> *vectorTarget2(0);
for(int iFace = 0; iFace < _nfaceQuad; iFace++) {
int indexLambda;
switch(iFace) {
case(0):
indexLambda = 5;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorV;
break;
case(1):
indexLambda = 3;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorW;
break;
case(2):
indexLambda = 1;
vectorTarget1 = &lkVectorV;
vectorTarget2 = &lkVectorW;
break;
case(3):
indexLambda = 0;
vectorTarget1 = &lkVectorV;
vectorTarget2 = &lkVectorW;
break;
case(4):
indexLambda = 2;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorW;
break;
case(5):
indexLambda = 4;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorV;
break;
}
for(int index1 = 0; index1 < _pOrderFace1[iFace] - 1; index1++) {
for(int index2 = 0; index2 < _pOrderFace2[iFace] - 1; index2++) {
faceBasis[indexFaceFunction] = lambda[indexLambda] *
(*vectorTarget1)[index1] *
(*vectorTarget2)[index2];
indexFaceFunction++;
}
}
}
// bubble shape functions:
int indexBubbleBasis = 0;
for(int ipb1 = 0; ipb1 < _pb1 - 1; ipb1++) {
for(int ipb2 = 0; ipb2 < _pb2 - 1; ipb2++) {
for(int ipb3 = 0; ipb3 < _pb3 - 1; ipb3++) {
bubbleBasis[indexBubbleBasis] =
lkVectorU[ipb1] * lkVectorV[ipb2] * lkVectorW[ipb3];
indexBubbleBasis++;
}
}
}
}
void HierarchicalBasisH1Brick::_someProductGrad(
double const &u, double const &v, double const &w,
std::vector<double> &product,
std::vector<std::vector<double> > &gradientProduct,
std::vector<double> &lambda,
std::vector<std::vector<double> > &gradientLambda)
{
lambda[0] = _affineCoordinate(1, u, v, w);
lambda[1] = _affineCoordinate(2, u, v, w);
lambda[2] = _affineCoordinate(3, u, v, w);
lambda[3] = _affineCoordinate(4, u, v, w);
lambda[4] = _affineCoordinate(5, u, v, w);
lambda[5] = _affineCoordinate(6, u, v, w);
gradientLambda[0][0] = 0.5;
gradientLambda[1][0] = -0.5;
gradientLambda[2][1] = 0.5;
gradientLambda[3][1] = -0.5;
gradientLambda[4][2] = 0.5;
gradientLambda[5][2] = -0.5;
product[0] = lambda[3] * lambda[5];
product[1] = lambda[1] * lambda[5];
product[2] = lambda[1] * lambda[3];
product[3] = lambda[0] * lambda[5];
product[4] = lambda[3] * lambda[0];
product[5] = lambda[2] * lambda[5];
product[6] = lambda[2] * lambda[0];
product[7] = lambda[2] * lambda[1];
product[8] = lambda[3] * lambda[4];
product[9] = lambda[4] * lambda[1];
product[10] = lambda[4] * lambda[0];
product[11] = lambda[4] * lambda[2];
gradientProduct[0][1] = -0.5 * lambda[5];
gradientProduct[0][2] = -0.5 * lambda[3];
gradientProduct[1][0] = -0.5 * lambda[5];
gradientProduct[1][2] = -0.5 * lambda[1];
gradientProduct[2][0] = -0.5 * lambda[3];
gradientProduct[2][1] = -0.5 * lambda[1];
gradientProduct[3][0] = 0.5 * lambda[5];
gradientProduct[3][2] = -0.5 * lambda[0];
gradientProduct[4][0] = 0.5 * lambda[3];
gradientProduct[4][1] = -0.5 * lambda[0];
gradientProduct[5][1] = 0.5 * lambda[5];
gradientProduct[5][2] = -0.5 * lambda[2];
gradientProduct[6][0] = 0.5 * lambda[2];
gradientProduct[6][1] = 0.5 * lambda[0];
gradientProduct[7][0] = -0.5 * lambda[2];
gradientProduct[7][1] = 0.5 * lambda[1];
gradientProduct[8][1] = -0.5 * lambda[4];
gradientProduct[8][2] = 0.5 * lambda[3];
gradientProduct[9][0] = -0.5 * lambda[4];
gradientProduct[9][2] = 0.5 * lambda[1];
gradientProduct[10][0] = 0.5 * lambda[4];
gradientProduct[10][2] = 0.5 * lambda[0];
gradientProduct[11][1] = 0.5 * lambda[4];
gradientProduct[11][2] = 0.5 * lambda[2];
}
void HierarchicalBasisH1Brick::generateGradientBasis(
double const &u, double const &v, double const &w,
std::vector<std::vector<double> > &gradientVertex,
std::vector<std::vector<double> > &gradientEdge,
std::vector<std::vector<double> > &gradientFace,
std::vector<std::vector<double> > &gradientBubble)
{
std::vector<double> product(12, 0);
std::vector<std::vector<double> > gradientProduct(12,
std::vector<double>(3, 0));
std::vector<double> lambda(6, 0);
std::vector<std::vector<double> > gradientLambda(6,
std::vector<double>(3, 0));
HierarchicalBasisH1Brick::_someProductGrad(u, v, w, product, gradientProduct,
lambda, gradientLambda);
// vertex gradient:
for(int it = 0; it < 3; it++) {
gradientVertex[0][it] =
gradientLambda[1][it] * product[0] + gradientProduct[0][it] * lambda[1];
gradientVertex[1][it] =
gradientLambda[0][it] * product[0] + gradientProduct[0][it] * lambda[0];
gradientVertex[2][it] =
gradientLambda[0][it] * product[5] + gradientProduct[5][it] * lambda[0];
gradientVertex[3][it] =
gradientLambda[1][it] * product[5] + gradientProduct[5][it] * lambda[1];
gradientVertex[4][it] =
gradientLambda[1][it] * product[8] + gradientProduct[8][it] * lambda[1];
gradientVertex[5][it] =
gradientLambda[0][it] * product[8] + gradientProduct[8][it] * lambda[0];
gradientVertex[6][it] =
gradientLambda[0][it] * product[11] + gradientProduct[11][it] * lambda[0];
gradientVertex[7][it] =
gradientLambda[1][it] * product[11] + gradientProduct[11][it] * lambda[1];
}
std::vector<double> lkVectorU(_pb1 - 1);
std::vector<double> lkVectorV(_pb2 - 1);
std::vector<double> lkVectorW(_pb3 - 1);
std::vector<std::vector<double> > dlkVectorU(_pb1 - 1,
std::vector<double>(3, 0.));
std::vector<std::vector<double> > dlkVectorV(_pb2 - 1,
std::vector<double>(3, 0.));
std::vector<std::vector<double> > dlkVectorW(_pb3 - 1,
std::vector<double>(3, 0.));
for(int it = 2; it <= _pb1; it++) {
lkVectorU[it - 2] = OrthogonalPoly::EvalLobatto(it, u);
dlkVectorU[it - 2][0] = OrthogonalPoly::EvalDLobatto(it, u);
}
for(int it = 2; it <= _pb2; it++) {
lkVectorV[it - 2] = OrthogonalPoly::EvalLobatto(it, v);
dlkVectorV[it - 2][1] = OrthogonalPoly::EvalDLobatto(it, v);
}
for(int it = 2; it <= _pb3; it++) {
lkVectorW[it - 2] = OrthogonalPoly::EvalLobatto(it, w);
dlkVectorW[it - 2][2] = OrthogonalPoly::EvalDLobatto(it, w);
}
// edge gradient:
int indexEdgeBasis = 0;
std::vector<double> *vectorTarget1(0);
std::vector<std::vector<double> > *dvectorTarget1(0);
for(int iEdge = 0; iEdge < _nedge; iEdge++) {
switch(iEdge) {
case(0):
case(5):
case(8):
case(11):
vectorTarget1 = &lkVectorU;
dvectorTarget1 = &dlkVectorU;
break;
case(1):
case(3):
case(9):
case(10):
vectorTarget1 = &lkVectorV;
dvectorTarget1 = &dlkVectorV;
break;
case(2):
case(4):
case(6):
case(7):
vectorTarget1 = &lkVectorW;
dvectorTarget1 = &dlkVectorW;
break;
}
for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] - 1;
indexEdgeFunc++) {
for(int it = 0; it < 3; it++) {
gradientEdge[indexEdgeBasis][it] =
(*dvectorTarget1)[indexEdgeFunc][it] * product[iEdge] +
(*vectorTarget1)[indexEdgeFunc] * gradientProduct[iEdge][it];
}
indexEdgeBasis++;
}
}
// face gradient:
int indexFaceFunction = 0;
std::vector<double> *vectorTarget2(0);
std::vector<std::vector<double> > *dvectorTarget2(0);
for(int iFace = 0; iFace < _nfaceQuad; iFace++) {
int indexLambda;
switch(iFace) {
case(0):
indexLambda = 5;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorV;
dvectorTarget1 = &dlkVectorU;
dvectorTarget2 = &dlkVectorV;
break;
case(1):
indexLambda = 3;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorW;
dvectorTarget1 = &dlkVectorU;
dvectorTarget2 = &dlkVectorW;
break;
case(2):
indexLambda = 1;
vectorTarget1 = &lkVectorV;
vectorTarget2 = &lkVectorW;
dvectorTarget1 = &dlkVectorV;
dvectorTarget2 = &dlkVectorW;
break;
case(3):
indexLambda = 0;
vectorTarget1 = &lkVectorV;
vectorTarget2 = &lkVectorW;
dvectorTarget1 = &dlkVectorV;
dvectorTarget2 = &dlkVectorW;
break;
case(4):
indexLambda = 2;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorW;
dvectorTarget1 = &dlkVectorU;
dvectorTarget2 = &dlkVectorW;
break;
case(5):
indexLambda = 4;
vectorTarget1 = &lkVectorU;
vectorTarget2 = &lkVectorV;
dvectorTarget1 = &dlkVectorU;
dvectorTarget2 = &dlkVectorV;
break;
}
for(int index1 = 0; index1 < _pOrderFace1[iFace] - 1; index1++) {
for(int index2 = 0; index2 < _pOrderFace2[iFace] - 1; index2++) {
for(int it = 0; it < 3; it++) {
gradientFace[indexFaceFunction][it] =
gradientLambda[indexLambda][it] * (*vectorTarget1)[index1] *
(*vectorTarget2)[index2] +
lambda[indexLambda] * (*dvectorTarget1)[index1][it] *
(*vectorTarget2)[index2] +
lambda[indexLambda] * (*vectorTarget1)[index1] *
(*dvectorTarget2)[index2][it];
}
indexFaceFunction++;
}
}
}
// bubble shape functions:
int indexBubbleBasis = 0;
for(int ipb1 = 0; ipb1 < _pb1 - 1; ipb1++) {
for(int ipb2 = 0; ipb2 < _pb2 - 1; ipb2++) {
for(int ipb3 = 0; ipb3 < _pb3 - 1; ipb3++) {
gradientBubble[indexBubbleBasis][0] =
dlkVectorU[ipb1][0] * lkVectorV[ipb2] * lkVectorW[ipb3];
gradientBubble[indexBubbleBasis][1] =
lkVectorU[ipb1] * dlkVectorV[ipb2][1] * lkVectorW[ipb3];
gradientBubble[indexBubbleBasis][2] =
lkVectorU[ipb1] * lkVectorV[ipb2] * dlkVectorW[ipb3][2];
indexBubbleBasis++;
}
}
}
}
void HierarchicalBasisH1Brick::orientEdge(int const &flagOrientation,
int const &edgeNumber,
std::vector<double> &edgeBasis)
{
if(flagOrientation == -1) {
int constant1 = 0;
int constant2 = 0;
for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
constant2 = constant2 - 1;
constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
for(int k = constant1; k <= constant2; k++) {
if((k - constant1) % 2 != 0) { edgeBasis[k] = edgeBasis[k] * (-1); }
}
}
}
void HierarchicalBasisH1Brick::orientEdge(
int const &flagOrientation, int const &edgeNumber,
std::vector<std::vector<double> > &gradientEdge)
{
if(flagOrientation == -1) {
int constant1 = 0;
int constant2 = 0;
for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
constant2 = constant2 - 1;
constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
for(int k = constant1; k <= constant2; k++) {
if((k - constant1) % 2 != 0) {
gradientEdge[k][0] = gradientEdge[k][0] * (-1);
gradientEdge[k][1] = gradientEdge[k][1] * (-1);
gradientEdge[k][2] = gradientEdge[k][2] * (-1);
}
}
}
}
void HierarchicalBasisH1Brick::orientFace(double const &u, double const &v,
double const &w, int const &flag1,
int const &flag2, int const &flag3,
int const &faceNumber,
std::vector<double> &faceBasis)
{
if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
int iterator = 0;
for(int i = 0; i < faceNumber; i++) {
iterator += (_pOrderFace1[i] - 1) * (_pOrderFace2[i] - 1);
}
if(flag3 == 1) {
for(int it1 = 2; it1 <= _pOrderFace1[faceNumber]; it1++) {
for(int it2 = 2; it2 <= _pOrderFace2[faceNumber]; it2++) {
int impactFlag1 = 1;
int impactFlag2 = 1;
if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
faceBasis[iterator] = faceBasis[iterator] * impactFlag1 * impactFlag2;
iterator++;
}
}
}
else {
double lambda = 0;
double var1 = 0;
double var2 = 0;
switch(faceNumber) {
case(0):
lambda = _affineCoordinate(6, u, v, w);
var1 = u;
var2 = v;
break;
case(1):
lambda = _affineCoordinate(4, u, v, w);
var1 = u;
var2 = w;
break;
case(2):
lambda = _affineCoordinate(2, u, v, w);
var1 = v;
var2 = w;
break;
case(3):
lambda = _affineCoordinate(1, u, v, w);
var1 = v;
var2 = w;
break;
case(4):
lambda = _affineCoordinate(3, u, v, w);
var1 = u;
var2 = w;
break;
case(5):
lambda = _affineCoordinate(5, u, v, w);
var1 = u;
var2 = v;
break;
}
std::vector<double> lkVector1(_pOrderFace1[faceNumber] - 1);
std::vector<double> lkVector2(_pOrderFace2[faceNumber] - 1);
for(int it = 2; it <= _pOrderFace1[faceNumber]; it++) {
lkVector1[it - 2] = OrthogonalPoly::EvalLobatto(it, var1);
}
for(int it = 2; it <= _pOrderFace2[faceNumber]; it++) {
lkVector2[it - 2] = OrthogonalPoly::EvalLobatto(it, var2);
}
for(int it1 = 2; it1 <= _pOrderFace2[faceNumber]; it1++) {
for(int it2 = 2; it2 <= _pOrderFace1[faceNumber]; it2++) {
int impactFlag1 = 1;
int impactFlag2 = 1;
if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
faceBasis[iterator] = lambda * lkVector1[it2 - 2] *
lkVector2[it1 - 2] * impactFlag1 * impactFlag2;
iterator++;
}
}
}
}
}
void HierarchicalBasisH1Brick::orientFace(
double const &u, double const &v, double const &w, int const &flag1,
int const &flag2, int const &flag3, int const &faceNumber,
std::vector<std::vector<double> > &gradientFace, std::string typeFunction)
{
if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
int iterator = 0;
for(int i = 0; i < faceNumber; i++) {
iterator += (_pOrderFace1[i] - 1) * (_pOrderFace2[i] - 1);
}
if(flag3 == 1) {
for(int it1 = 2; it1 <= _pOrderFace1[faceNumber]; it1++) {
for(int it2 = 2; it2 <= _pOrderFace2[faceNumber]; it2++) {
int impactFlag1 = 1;
int impactFlag2 = 1;
if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
gradientFace[iterator][0] =
gradientFace[iterator][0] * impactFlag1 * impactFlag2;
gradientFace[iterator][1] =
gradientFace[iterator][1] * impactFlag1 * impactFlag2;
gradientFace[iterator][2] =
gradientFace[iterator][2] * impactFlag1 * impactFlag2;
iterator++;
}
}
}
else {
std::vector<double> uvw(3);
uvw[0] = u;
uvw[1] = v;
uvw[2] = w;
double lambda = 0;
int var1 = 0;
int var2 = 0;
std::vector<double> gradientLambda(3, 0);
switch(faceNumber) {
case(0):
lambda = _affineCoordinate(6, u, v, w);
var1 = 0;
var2 = 1;
gradientLambda[2] = -0.5;
break;
case(1):
lambda = _affineCoordinate(4, u, v, w);
var1 = 0;
var2 = 2;
gradientLambda[1] = -0.5;
break;
case(2):
lambda = _affineCoordinate(2, u, v, w);
var1 = 1;
var2 = 2;
gradientLambda[0] = -0.5;
break;
case(3):
lambda = _affineCoordinate(1, u, v, w);
var1 = 1;
var2 = 2;
gradientLambda[0] = 0.5;
break;
case(4):
lambda = _affineCoordinate(3, u, v, w);
var1 = 0;
var2 = 2;
gradientLambda[1] = 0.5;
break;
case(5):
lambda = _affineCoordinate(5, u, v, w);
var1 = 0;
var2 = 1;
gradientLambda[2] = 0.5;
break;
}
std::vector<double> lkVector1(_pOrderFace1[faceNumber] - 1);
std::vector<double> lkVector2(_pOrderFace2[faceNumber] - 1);
std::vector<std::vector<double> > dlkVector1(_pOrderFace1[faceNumber] - 1,
std::vector<double>(3, 0.));
std::vector<std::vector<double> > dlkVector2(_pOrderFace2[faceNumber] - 1,
std::vector<double>(3, 0.));
for(int it = 2; it <= _pOrderFace1[faceNumber]; it++) {
lkVector1[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[var1]);
dlkVector1[it - 2][var1] = OrthogonalPoly::EvalDLobatto(it, uvw[var1]);
}
for(int it = 2; it <= _pOrderFace2[faceNumber]; it++) {
lkVector2[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[var2]);
dlkVector2[it - 2][var2] = OrthogonalPoly::EvalDLobatto(it, uvw[var2]);
}
for(int it1 = 2; it1 <= _pOrderFace2[faceNumber]; it1++) {
for(int it2 = 2; it2 <= _pOrderFace1[faceNumber]; it2++) {
int impactFlag1 = 1;
int impactFlag2 = 1;
if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
for(int itVector = 0; itVector < 3; itVector++) {
gradientFace[iterator][itVector] =
(gradientLambda[itVector] * lkVector1[it2 - 2] *
lkVector2[it1 - 2] +
lambda * dlkVector1[it2 - 2][itVector] * lkVector2[it1 - 2] +
lambda * lkVector1[it2 - 2] * dlkVector2[it1 - 2][itVector]) *
impactFlag1 * impactFlag2;
}
iterator++;
}
}
}
}
}