/*********************************************************************/ /* File: highorderfe.cc */ /* Author: Joachim Schoeberl */ /* Date: 28. Oct. 2000 */ /*********************************************************************/ /* High Order Finite Elements */ #include namespace ngfem { using namespace ngfem; FE_SegmP :: FE_SegmP (int ap, bool ahierarchical) : NodalFiniteElement (1, ET_SEGM, ap+1, ap) { p = ap; hierarchical = ahierarchical; if (!ipdata.Size()) CalcIPData(ET_SEGM, ipdata); } FE_SegmP :: ~FE_SegmP() { /* for (int i = 0; i < ipdata.Size(); i++) delete ipdata[i]; */ } void FE_SegmP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { double x[3]; x[0] = ip(0); x[1] = 1 - x[0]; int i = 0; int j, k, l, m, n; int lami[2]; if (hierarchical) { // vertex shape functions shape(0) = x[0]; shape(p) = x[1]; i = 1; for (lami[0] = 1; lami[0] < p; lami[0]++) { lami[1] = p - lami[0]; double fac = 1; for (n = 0; n < 2; n++) for (m = 0; m < lami[n]; m++) fac *= (p * x[n] - m) / double (lami[n] - m); shape(i) = fac; i++; } } else { // all nodal i = 0; for (lami[0] = 0; lami[0] <= p; lami[0]++) { lami[1] = p - lami[0]; double fac = 1; for (n = 0; n < 2; n++) for (m = 0; m < lami[n]; m++) fac *= (p * x[n] - m) / double (lami[n] - m); shape(i) = fac; i++; } } } // ARRAY FE_TrigP::ipdata; FE_TrigP :: FE_TrigP (int ap) : NodalFiniteElement (2, ET_TRIG, (ap*ap + 3*ap+2) / 2, ap) { p = ap; hierarchical = 1; if (!ipdata.Size()) CalcIPData(ET_TRIG, ipdata); } FE_TrigP :: ~FE_TrigP() { ; } void FE_TrigP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { double x[3], px[3]; x[0] = ip(0); x[1] = ip(1); x[2] = 1 - x[0] - x[1]; int i, j, k, l, m, n; int lami[3]; for (i = 0; i < 3; i++) px[i] = p * x[i]; i = 0; for (lami[0] = 0; lami[0] <= p; lami[0]++) for (lami[1] = 0; lami[1] <= p-lami[0]; lami[1]++) { lami[2] = p - lami[0] - lami[1]; // check for linear shape functions int corner = -1; for (n = 0; n < 3; n++) if (lami[n] == p) corner = n; if (corner >= 0) { shape(i) = x[corner]; i++; continue; } double reducex = 0; for (n = 0; n < 3; n++) if (lami[n] == 0) reducex += x[n]; reducex -= 1; double fac = 1; for (n = 0; n < 3; n++) for (m = 0; m < lami[n]; m++) fac *= (px[n] + m * reducex) / double (m+1); shape(i) = fac; i++; } } #ifdef NONE // ARRAY FE_QuadP::ipdata; FE_QuadP :: FE_QuadP (int ap, int ap2) : segm(ap) { p = ap; if (ap2 == -1) p2 = p; else p2 = ap2; if (p == p2) segm2 = &segm; else segm2 = new FE_SegmP(p2); if (!ipdata.Size()) CalcIPData(ET_QUAD, ipdata); } FE_QuadP :: ~FE_QuadP() { int i; /* for (i = 1; i <= ipdata.Size(); i++) delete ipdata.Elem(i); */ if (segm2 != &segm) delete segm2; } int FE_QuadP :: GetNDof () const { return (p+1)*(p2+1); } void FE_QuadP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { // shape.SetSize(GetNDof()); IntegrationPoint ip1(ip.Point()[0], 0, 0, 1); IntegrationPoint ip2(ip.Point()[1], 0, 0, 1); Vector<> shape1 (segm.GetNDof()); Vector<> shape2 (segm2->GetNDof()); segm.CalcShape (ip1, shape1); segm2->CalcShape (ip2, shape2); int i = 0; int j, k; for (j = 1; j <= p2+1; j++) for (k = 1; k <= p+1; k++) { i++; shape.Elem(i) = shape1.Get(k) * shape2.Get(j); } } #endif FE_TetP :: FE_TetP (int ap) : NodalFiniteElement (3, ET_TET, (ap*ap*ap + 6 * ap * ap + 11 * ap + 6) / 6, ap) { p = ap; hierarchical = 1; if (!ipdata.Size()) CalcIPData(ET_TET, ipdata); } FE_TetP :: ~FE_TetP() { ; } /* void FE_TetP :: CalcShape (const IntegrationPoint & ip, Vector & shape, int comp) const { shape.SetSize(GetNDof()); double x[4]; x[0] = ip.Point()[0]; x[1] = ip.Point()[1]; x[2] = ip.Point()[2]; x[3] = 1 - x[0] - x[1] - x[2]; int i = 0; int j, k, l, m, n; int lami[4]; double xi[4]; for (lami[0] = 0; lami[0] <= p; lami[0]++) for (lami[1] = 0; lami[1] <= p-lami[0]; lami[1]++) for (lami[2] = 0; lami[2] <= p-lami[0]-lami[1]; lami[2]++) { lami[3] = p - lami[0] - lami[1] - lami[2]; for (n = 0; n < 4; n++) xi[n] = double(lami[n]) / p; double fac = 1; for (n = 0; n < 4; n++) for (m = 0; m < lami[n]; m++) fac *= (x[n] - double(m) / p) / (xi[n] - double(m) / p); i++; shape.Elem(i) = fac; } } */ void FE_TetP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { // shape.SetSize(GetNDof()); double x[4], px[4]; x[0] = ip.Point()[0]; x[1] = ip.Point()[1]; x[2] = ip.Point()[2]; x[3] = 1 - x[0] - x[1] - x[2]; int i, j, k, l, m, n; int lami[4]; for (i = 0; i < 4; i++) px[i] = p * x[i]; i = 0; for (lami[0] = 0; lami[0] <= p; lami[0]++) for (lami[1] = 0; lami[1] <= p-lami[0]; lami[1]++) for (lami[2] = 0; lami[2] <= p-lami[0]-lami[1]; lami[2]++) { lami[3] = p - lami[0] - lami[1] - lami[2]; int corner = -1; for (n = 0; n < 4; n++) if (lami[n] == p) corner = n; if (corner >= 0) { shape(i) = x[corner]; i++; continue; } double redx = 0; for (n = 0; n < 4; n++) if (lami[n] == 0) redx += x[n]; redx -= 1; double fac = 1; for (n = 0; n < 4; n++) for (m = 0; m < lami[n]; m++) fac *= (px[n] + redx * m) / double(m+1); shape(i) = fac; i++; } /* (*testout) << "p = " << ip.Point()[0] << ", " << ip.Point()[1] << ", " << ip.Point()[2] << ": " << endl << shape << endl; */ } #ifdef NONE // ARRAY FE_PrismP::ipdata; FE_PrismP :: FE_PrismP (int ap) { p = ap; if (!ipdata.Size()) CalcIPData(ET_PRISM, ipdata); } FE_PrismP :: ~FE_PrismP() { ; } int FE_PrismP :: GetNDof () const { return (p+1) * (p*p + 3*p+2) / 2; } void FE_PrismP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { // shape.SetSize(GetNDof()); FE_TrigP trig(p); FE_SegmP segm(p); IntegrationPoint ip1(ip.Point()[0], ip.Point()[1], 0, 1); IntegrationPoint ip2(ip.Point()[2], 0, 0, 1); Vector<> shape1 (trig.GetNDof()); Vector<> shape2 (segm.GetNDof()); trig.CalcShape (ip1, shape1); segm.CalcShape (ip2, shape2); int i = 0; int j, k; for (j = 1; j <= shape2.Size(); j++) for (k = 1; k <= shape1.Size(); k++) { shape(i) = shape1.Get(k) * shape2.Get(j); i++; } } // ARRAY FE_HexP::ipdata; FE_HexP :: FE_HexP (int ap) { p = ap; if (!ipdata.Size()) CalcIPData(ET_HEX, ipdata); } FE_HexP :: ~FE_HexP() { ; } int FE_HexP :: GetNDof () const { return (p+1) * (p+1) * (p+1); } void FE_HexP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { // shape.SetSize(GetNDof()); FE_SegmP segm(p); IntegrationPoint ip1(ip.Point()[0], 0, 0, 1); IntegrationPoint ip2(ip.Point()[1], 0, 0, 1); IntegrationPoint ip3(ip.Point()[2], 0, 0, 1); Vector<> shape1 (segm.GetNDof()); Vector<> shape2 (segm.GetNDof()); Vector<> shape3 (segm.GetNDof()); segm.CalcShape (ip1, shape1); segm.CalcShape (ip2, shape2); segm.CalcShape (ip3, shape3); int i = 0; int j, k, l; for (j = 1; j <= shape3.Size(); j++) for (k = 1; k <= shape2.Size(); k++) for (l = 1; l <= shape1.Size(); l++) { shape(i) = shape1.Get(l) * shape2.Get(k) * shape3.Get(j); i++; } } #endif FE_Augmented_SegmP :: FE_Augmented_SegmP (int ap) : NodalFiniteElement (1, ET_SEGM, ap+1+2*(ap-1), ap) { p = ap; if (!ipdata.Size()) { CalcIPData(ET_SEGM, ipdata); // compute unification matrix int i, j, k; int ndreg = p+1; int nd = ndreg + 2*(p-1); Matrix<> mat1(ndreg, ndreg); Matrix<> imat1(ndreg, ndreg); Matrix<> mat2(ndreg, nd-ndreg); Matrix<> mat_unify(ndreg, nd-ndreg); Vector<> shape(nd); int ii = 0; for (i = 0; i <= p; i++) { IntegrationPoint ip(double(i)/p, 0, 0, 0); CalcShape (ip, shape); for (k = 0; k < ndreg; k++) mat1(ii, k) = shape(k); for (k = 0; k < nd-ndreg; k++) mat2(ii, k) = shape(k+ndreg); ii++; } CalcInverse (mat1, imat1); mat_unify = imat1 * mat2; (*testout) << "mat1 = " << endl << mat1 << endl << "mat2 = " << endl << mat2 << endl << "unify = " << endl << mat_unify << endl; } } FE_Augmented_SegmP :: ~FE_Augmented_SegmP() { /* for (int i = 0; i < ipdata.Size(); i++) delete ipdata[i]; */ } void FE_Augmented_SegmP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { double x[3]; x[0] = ip(0); x[1] = 1 - x[0]; int i = 0; int j, k, l, m, n; int lami[2]; // vertex shape functions shape(0) = x[0]; shape(p) = x[1]; i = 1; for (lami[0] = 1; lami[0] < p; lami[0]++) { lami[1] = p - lami[0]; double fac = 1; for (n = 0; n < 2; n++) for (m = 0; m < lami[n]; m++) fac *= (p * x[n] - m) / double (lami[n] - m); shape(i) = fac; i++; } i = p+1; for (j = 0; j < 2; j++) { for (k = 2; k <= p; k++) { shape(i) = pow (x[j], k); i++; } } } // ARRAY FE_Augmented_TrigP::ipdata; FE_Augmented_TrigP :: FE_Augmented_TrigP (int ap) : NodalFiniteElement (2, ET_TRIG, (ap*ap + 3*ap+2) / 2 + 3*(ap-1), ap), mat_unify ( (ap*ap + 3*ap+2) / 2, 3*(ap-1)) { cout << " allocate augmented TrigP" << endl; p = ap; if (!ipdata.Size()) { CalcIPData(ET_TRIG, ipdata); // compute unification matrix int i, j, k; int ndreg = (p*p + 3*p+2) / 2; int nd = ndreg + 3*(p-1); Matrix<> mat1(ndreg, ndreg); Matrix<> imat1(ndreg, ndreg); Matrix<> mat2(ndreg, nd-ndreg); Vector<> shape(nd); int ii = 0; for (i = 0; i <= p; i++) for (j = 0; j <= p-i; j++) { IntegrationPoint ip(double(i)/p, double(j)/p, 0, 0); CalcShape (ip, shape); for (k = 0; k < ndreg; k++) mat1(ii, k) = shape(k); for (k = 0; k < nd-ndreg; k++) mat2(ii, k) = shape(k+ndreg); ii++; } CalcInverse (mat1, imat1); mat_unify = imat1 * mat2; (*testout) << "mat1 = " << endl << mat1 << endl << "mat2 = " << endl << mat2 << endl << "unify = " << endl << mat_unify << endl; } } FE_Augmented_TrigP :: ~FE_Augmented_TrigP() { ; } void FE_Augmented_TrigP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { double x[3], px[3]; x[0] = ip(0); x[1] = ip(1); x[2] = 1 - x[0] - x[1]; int i, j, k, l, m, n; int lami[3]; for (i = 0; i < 3; i++) px[i] = p * x[i]; i = 0; for (lami[0] = 0; lami[0] <= p; lami[0]++) for (lami[1] = 0; lami[1] <= p-lami[0]; lami[1]++) { lami[2] = p - lami[0] - lami[1]; // check for linear shape functions int corner = -1; for (n = 0; n < 3; n++) if (lami[n] == p) corner = n; if (corner >= 0) { shape(i) = x[corner]; i++; continue; } double reducex = 0; for (n = 0; n < 3; n++) if (lami[n] == 0) reducex += x[n]; reducex -= 1; double fac = 1; for (n = 0; n < 3; n++) for (m = 0; m < lami[n]; m++) fac *= (px[n] + m * reducex) / double (m+1); shape(i) = fac; i++; } for (j = 0; j < 3; j++) { for (k = 2; k <= p; k++) { shape(i) = pow (x[j], k); i++; } } } void FE_Augmented_TrigP :: Unify (FlatVector<> & uelem) const { int ndreg = (p*p + 3*p+2) / 2; int nd = ndreg + 3*(p-1); FlatVector<> uelem1 (ndreg, &uelem(0)); FlatVector<> uelem2 (nd-ndreg, &uelem(ndreg)); uelem1 += mat_unify * uelem2; uelem2 = 0.0; } void FE_Augmented_TrigP :: UnifyTrans (FlatVector<> & uelem) const { int ndreg = (p*p + 3*p+2) / 2; int nd = ndreg + 3*(p-1); FlatVector<> uelem1 (ndreg, &uelem(0)); FlatVector<> uelem2 (nd-ndreg, &uelem(ndreg)); uelem2 = Trans (mat_unify) * uelem1; } FE_Augmented_TetP :: FE_Augmented_TetP (int ap) : NodalFiniteElement (3, ET_TET, (ap*ap*ap + 6 * ap * ap + 11 * ap + 6) / 6 + 4 * (ap-1), ap), mat_unify ( (ap*ap*ap + 6 * ap * ap + 11 * ap + 6) / 6, 4*(ap-1)) { p = ap; if (!ipdata.Size()) { CalcIPData(ET_TET, ipdata); // compute unification matrix int i, j, k, l; int ndreg = (p*p*p + 6 * p * p + 11 * p + 6) / 6; int nd = ndreg + 4*(p-1); Matrix<> mat1(ndreg, ndreg); Matrix<> imat1(ndreg, ndreg); Matrix<> mat2(ndreg, nd-ndreg); Vector<> shape(nd); int ii = 0; for (i = 0; i <= p; i++) for (j = 0; j <= p-i; j++) for (l = 0; l <= p-i-j; l++) { IntegrationPoint ip(double(i)/p, double(j)/p, double(l)/p, 0); CalcShape (ip, shape); for (k = 0; k < ndreg; k++) mat1(ii, k) = shape(k); for (k = 0; k < nd-ndreg; k++) mat2(ii, k) = shape(k+ndreg); ii++; } CalcInverse (mat1, imat1); mat_unify = imat1 * mat2; (*testout) << "mat1 = " << endl << mat1 << endl << "mat2 = " << endl << mat2 << endl << "unify = " << endl << mat_unify << endl; } } FE_Augmented_TetP :: ~FE_Augmented_TetP() { ; } void FE_Augmented_TetP :: CalcShape (const IntegrationPoint & ip, FlatVector<> shape) const { // shape.SetSize(GetNDof()); double x[4], px[4]; x[0] = ip(0); x[1] = ip(1); x[2] = ip(2); x[3] = 1 - x[0] - x[1] - x[2]; int i, j, k, l, m, n; int lami[4]; for (i = 0; i < 4; i++) px[i] = p * x[i]; i = 0; for (lami[0] = 0; lami[0] <= p; lami[0]++) for (lami[1] = 0; lami[1] <= p-lami[0]; lami[1]++) for (lami[2] = 0; lami[2] <= p-lami[0]-lami[1]; lami[2]++) { lami[3] = p - lami[0] - lami[1] - lami[2]; int corner = -1; for (n = 0; n < 4; n++) if (lami[n] == p) corner = n; if (corner >= 0) { shape(i) = x[corner]; i++; continue; } double redx = 0; for (n = 0; n < 4; n++) if (lami[n] == 0) redx += x[n]; redx -= 1; double fac = 1; for (n = 0; n < 4; n++) for (m = 0; m < lami[n]; m++) fac *= (px[n] + redx * m) / double(m+1); shape(i) = fac; i++; } for (j = 0; j < 4; j++) { for (k = 2; k <= p; k++) { shape(i) = pow (x[j], k); i++; } } } void FE_Augmented_TetP :: Unify (FlatVector<> & uelem) const { int ndreg = (p*p*p + 6 * p * p + 11 * p + 6) / 6; int nd = ndreg + 4*(p-1); FlatVector<> uelem1 (ndreg, &uelem(0)); FlatVector<> uelem2 (nd-ndreg, &uelem(ndreg)); uelem1 += mat_unify * uelem2; uelem2 = 0.0; } void FE_Augmented_TetP :: UnifyTrans (FlatVector<> & uelem) const { int ndreg = (p*p*p + 6 * p * p + 11 * p + 6) / 6; int nd = ndreg + 4*(p-1); FlatVector<> uelem1 (ndreg, &uelem(0)); FlatVector<> uelem2 (nd-ndreg, &uelem(ndreg)); uelem2 = Trans (mat_unify) * uelem1; } }