// Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle // // See the LICENSE.txt file for license information. Please report all // bugs and problems to the public mailing list <gmsh@onelab.info>. #include "elasticityTerm.h" #include "Numeric.h" // The SElement (Solver element) that has been sent to the function // contains 2 enrichments, that can enrich both shape and test functions void elasticityTerm::createData(MElement *e) const { if (_data.find(e->getTypeForMSH()) != _data.end())return; elasticityDataAtGaussPoint d; int nbSF = e->getNumShapeFunctions(); int integrationOrder = 2 * (e->getPolynomialOrder() - 1) ; int npts; IntPt *GP; double gs[100][3]; e->getIntegrationPoints(integrationOrder, &npts, &GP); for (int i=0;i<npts;i++){ fullMatrix<double> a(nbSF,3); const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2]; const double weight = GP[i].weight; e->getGradShapeFunctions(u,v,w,gs); for (int j=0;j < nbSF;j++){ a(j,0) = gs[j][0]; a(j,1) = gs[j][1]; a(j,2) = gs[j][2]; } d.gradSF.push_back(a); d.u.push_back(u); d.v.push_back(v); d.w.push_back(w); d.weight.push_back(weight); } // printf("coucou creation of a data for %d with %d points\n",e->getTypeForMSH(),npts); _data[e->getTypeForMSH()] = d; } void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const { MElement *e = se->getMeshElement(); createData(e); int nbSF = e->getNumShapeFunctions(); elasticityDataAtGaussPoint &d = _data[e->getTypeForMSH()]; int npts = d.u.size(); m.setAll(0.); double FACT = _E / (1 + _nu); double C11 = FACT * (1 - _nu) / (1 - 2 * _nu); double C12 = FACT * _nu / (1 - 2 * _nu); double C44 = (C11 - C12) / 2; const double C[6][6] = { {C11, C12, C12, 0, 0, 0}, {C12, C11, C12, 0, 0, 0}, {C12, C12, C11, 0, 0, 0}, { 0, 0, 0, C44, 0, 0}, { 0, 0, 0, 0, C44, 0}, { 0, 0, 0, 0, 0, C44} }; fullMatrix<double> H(6, 6); fullMatrix<double> B(6, 3 * nbSF); fullMatrix<double> BTH(3 * nbSF, 6); fullMatrix<double> BT(3 * nbSF, 6); for (int i = 0; i < 6; i++) for (int j = 0; j < 6; j++) H(i, j) = C[i][j]; double jac[3][3],invjac[3][3],Grads[100][3]; for (int i = 0; i < npts; i++){ const double weight = d.weight[i]; const fullMatrix<double> &grads = d.gradSF[i]; const double detJ = e->getJacobian(grads, jac); inv3x3(jac, invjac); for (int j = 0; j < nbSF; j++){ Grads[j][0] = invjac[0][0] * grads(j,0) + invjac[0][1] * grads(j,1) + invjac[0][2] * grads(j,2); Grads[j][1] = invjac[1][0] * grads(j,0) + invjac[1][1] * grads(j,1) + invjac[1][2] * grads(j,2); Grads[j][2] = invjac[2][0] * grads(j,0) + invjac[2][1] * grads(j,1) + invjac[2][2] * grads(j,2); } B.setAll(0.); BT.setAll(0.); if (se->getShapeEnrichement() == se->getTestEnrichement()){ // printf("coucou\n"); for (int j = 0; j < nbSF; j++){ BT(j, 0) = B(0, j) = Grads[j][0]; BT(j, 3) = B(3, j) = Grads[j][1]; BT(j, 5) = B(5, j) = Grads[j][2]; BT(j + nbSF, 1) = B(1, j + nbSF) = Grads[j][1]; BT(j + nbSF, 3) = B(3, j + nbSF) = Grads[j][0]; BT(j + nbSF, 4) = B(4, j + nbSF) = Grads[j][2]; BT(j + 2 * nbSF, 2) = B(2, j + 2 * nbSF) = Grads[j][2]; BT(j + 2 * nbSF, 4) = B(4, j + 2 * nbSF) = Grads[j][1]; BT(j + 2 * nbSF, 5) = B(5, j + 2 * nbSF) = Grads[j][0]; } } else{ /* se->gradNodalTestFunctions (u, v, w, invjac,GradsT); for (int j = 0; j < nbSF; j++){ BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0]; BT(j, 3) = Grads[j][1]; B(3, j) = GradsT[j][1]; BT(j, 5) = Grads[j][2]; B(5, j) = GradsT[j][2]; BT(j + nbSF, 1) = Grads[j][1]; B(1, j + nbSF) = GradsT[j][1]; BT(j + nbSF, 3) = Grads[j][0]; B(3, j + nbSF) = GradsT[j][0]; BT(j + nbSF, 4) = Grads[j][2]; B(4, j + nbSF) = GradsT[j][2]; BT(j + 2 * nbSF, 2) = Grads[j][2]; B(2, j + 2 * nbSF) = GradsT[j][2]; BT(j + 2 * nbSF, 4) = Grads[j][1]; B(4, j + 2 * nbSF) = GradsT[j][1]; BT(j + 2 * nbSF, 5) = Grads[j][0]; B(5, j + 2 * nbSF) = GradsT[j][0]; } */ } BTH.setAll(0.); BTH.gemm(BT, H); m.gemm(BTH, B, weight * detJ, 1.); } return; } void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const { MElement *e = se->getMeshElement(); int nbSF = e->getNumShapeFunctions(); int integrationOrder = 2 * e->getPolynomialOrder(); int npts; IntPt *GP; double jac[3][3]; double ff[256]; e->getIntegrationPoints(integrationOrder, &npts, &GP); m.scale(0.); for (int i = 0; i < npts; i++){ const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2]; const double weight = GP[i].weight; const double detJ = e->getJacobian(u, v, w, jac); se->nodalTestFunctions(u, v, w, ff); for (int j = 0; j < nbSF; j++){ m(j) += _volumeForce.x() * ff[j] * weight * detJ * .5; m(j + nbSF) += _volumeForce.y() * ff[j] * weight * detJ * .5; m(j + 2 * nbSF) += _volumeForce.z() * ff[j] * weight * detJ * .5; } } } /* B_d is the deviatoric part of the FE strains K_{uu} = \int_\Omega B^T_d H B_d dv */ void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const { setPolynomialBasis(se); MElement *e = se->getMeshElement(); int integrationOrder = 4 * (e->getPolynomialOrder()) + 2; int npts; IntPt *GP; double jac[3][3]; double invjac[3][3]; double Grads[100][3]; double SF[100]; e->getIntegrationPoints(integrationOrder, &npts, &GP); m.setAll(0.); fullMatrix<double> KUU(3 * _sizeN,3 * _sizeN); fullMatrix<double> KUP(3 * _sizeN,_sizeM); fullMatrix<double> KUG(3 * _sizeN,_sizeM); fullMatrix<double> KPG(_sizeM,_sizeM); fullMatrix<double> KGG(_sizeM,_sizeM); fullMatrix<double> KPP(_sizeM,_sizeM); // stabilization double FACT = _E / ((1 + _nu)*(1.-2*_nu)); double C11 = FACT * (1 - _nu) ; double C12 = FACT * _nu ; double C44 = FACT* (1 - 2.*_nu)*.5; //double _K = 3*_E / (1 - 2 * _nu); // FIXME : PLANE STRESS !!! // FACT = _E / (1-_nu*_nu); // C11 = FACT; // C12 = _nu * FACT; // C44 = (1.-_nu)*.5*FACT; const double C[6][6] = { {C11, C12, C12, 0, 0, 0}, {C12, C11, C12, 0, 0, 0}, {C12, C12, C11, 0, 0, 0}, { 0, 0, 0, C44, 0, 0}, { 0, 0, 0, 0, C44, 0}, { 0, 0, 0, 0, 0, C44} }; fullMatrix<double> H (6, 6); fullMatrix<double> B (6, 3 * _sizeN); fullMatrix<double> BDEV (6, 3 * _sizeN); fullMatrix<double> BDIL (1, 3 * _sizeN); fullMatrix<double> BDILT (3 * _sizeN,1); fullMatrix<double> BTH(3 * _sizeN, 6); fullMatrix<double> BT(3 * _sizeN, 6); fullMatrix<double> trBTH(3 * _sizeN, 1); fullMatrix<double> MT(_sizeM, 1); fullMatrix<double> M(1,_sizeM); for (int i = 0; i < 6; i++) for (int j = 0; j < 6; j++) H(i, j) = C[i][j]; double _dimension = 2.0; for (int i = 0; i < npts; i++){ const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2]; const double weight = GP[i].weight; const double detJ = e->getJacobian(u, v, w, jac); inv3x3(jac, invjac); se->gradNodalShapeFunctions(u, v, w, invjac,Grads); e->getShapeFunctions (u,v,w,SF,_polyOrderM); B.setAll(0.); BT.setAll(0.); const double K = .0000002*weight * detJ * pow (detJ,2/_dimension); // the second detJ is for stabilization for (int j = 0; j < _sizeM; j++){ for (int k = 0; k < _sizeM; k++){ KPP(j, k) += (Grads[j][0] * Grads[k][0] + Grads[j][1] * Grads[k][1] + Grads[j][2] * Grads[k][2])* K; } } const double K2 = weight * detJ; for (int j = 0; j < _sizeN; j++){ for (int k = 0; k < _sizeM; k++){ KUP(j+0*_sizeN, k) += Grads[j][0] * SF[k] * K2; KUP(j+1*_sizeN, k) += Grads[j][1] * SF[k] * K2; KUP(j+2*_sizeN, k) += Grads[j][2] * SF[k] * K2; } } /* const double K3 = weight * detJ * _E/(2*(1+_nu)); for (int j = 0; j < _sizeN; j++){ for (int k = 0; k < _sizeN; k++){ const double fa = (Grads[j][0] * Grads[k][0] + Grads[j][1] * Grads[k][1] + Grads[j][2] * Grads[k][2]) * K3; KUU(j+0*_sizeN, k+0*_sizeN) += fa; KUU(j+1*_sizeN, k+1*_sizeN) += fa; KUU(j+2*_sizeN, k+2*_sizeN) += fa; } } */ for (int j = 0; j < _sizeN; j++){ // printf("%g %g %g\n",Grads[j][0],Grads[j][1],Grads[j][2]); BDILT(j, 0) = BDIL(0, j) = Grads[j][0]/_dimension; BT(j, 0) = B(0, j) = Grads[j][0];//-Grads[j][0]/_dimension; // BT(j, 1) = B(1, j) = -Grads[j][1]/_dimension; // BT(j, 2) = B(2, j) = -Grads[j][2]/_dimension; BT(j, 3) = B(3, j) = Grads[j][1]; BT(j, 4) = B(4, j) = Grads[j][2]; BDILT(j + _sizeN, 0) = BDIL(0, j + _sizeN) = Grads[j][1]/_dimension; // BT(j + _sizeN, 0) = B(0, j + _sizeN) = -Grads[j][0]/_dimension; BT(j + _sizeN, 1) = B(1, j + _sizeN) = Grads[j][1];//-Grads[j][1]/_dimension; // BT(j + _sizeN, 2) = B(2, j + _sizeN) = -Grads[j][2]/_dimension; BT(j + _sizeN, 3) = B(3, j + _sizeN) = Grads[j][0]; BT(j + _sizeN, 5) = B(5, j + _sizeN) = Grads[j][2]; BDILT(j + 2 * _sizeN, 0) = BDIL(0, j + 2 * _sizeN) = Grads[j][2]/_dimension; // BT(j + 2 * _sizeN, 0) = B(0, j + 2 * _sizeN) = -Grads[j][0]/_dimension; // BT(j + 2 * _sizeN, 1) = B(1, j + 2 * _sizeN) = -Grads[j][1]/_dimension; BT(j + 2 * _sizeN, 2) = B(2, j + 2 * _sizeN) = Grads[j][2];//-Grads[j][2]/_dimension; BT(j + 2 * _sizeN, 4) = B(4, j + 2 * _sizeN) = Grads[j][1]; BT(j + 2 * _sizeN, 5) = B(5, j + 2 * _sizeN) = Grads[j][0]; } for (int j = 0; j < _sizeM; j++){ MT(j,0) = M(0,j) = SF[j]; } // printf("BT (%d %d) x H (%d,%d) = BTH(%d,%d)\n",BT.size1(),BT.size2(), H.size1(),H.size2(),BTH.size1(),BTH.size2()); BTH.setAll(0.); BTH.gemm(BT, H); for (int j = 0; j < 3*_sizeN; j++){ trBTH(j,0) = BTH(j,0) + BTH(j,1) + BTH(j,2) ; } KUU.gemm(BTH, B, weight * detJ); // KUP.gemm(BDILT, M, weight * detJ); // KPG.gemm(MT, M, -weight * detJ); // KUG.gemm(trBTH, M, weight * detJ/_dimension); // KGG.gemm(MT, M, weight * detJ * (_K)/(_dimension*_dimension)); } /* 3*_sizeN _sizeM _sizeM KUU KUP KUG m = KUP^t KPP KPG KUG^T KPG^t KGG */ for (int i=0;i<3*_sizeN;i++){ for (int j=0;j<3*_sizeN;j++) m(i,j) = KUU(i,j); // assemble KUU for (int j=0;j<_sizeM;j++){ m(i,j+3*_sizeN) = KUP(i,j); // assemble KUP m(j+3*_sizeN,i) = KUP(i,j); // assemble KUP^T } for (int j=0;j<_sizeM;j++){ m(i,j+3*_sizeN+_sizeM) = KUG(i,j); // assemble KUG m(j+3*_sizeN+_sizeM,i) = KUG(i,j); // assemble KUG^t } } for (int i=0;i<_sizeM;i++){ for (int j=0;j<_sizeM;j++){ m(i+3*_sizeN,j+3*_sizeN+_sizeM) = KPG(i,j); // assemble KPG m(j+3*_sizeN+_sizeM,i+3*_sizeN) = KPG(i,j); // assemble KPG^t m(i+3*_sizeN+_sizeM,j+3*_sizeN+_sizeM) = KGG(i,j); // assemble KGG m(i+3*_sizeN,j+3*_sizeN) = KPP(i,j); // assemble KPP } } // m.print("Mixed"); return; }