From 9c20c79ae97bac353d9850f96f54733637a7d93a Mon Sep 17 00:00:00 2001 From: kmo Date: Thu, 20 Oct 2011 16:56:12 +0000 Subject: [PATCH] ..and finally the Qq/Qq-1 formulations git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1273 e10b68d5-8a6e-419e-a041-bce267b0401d --- .../NonlinearElasticityULMixed.C | 24 +++++++++++++++---- .../NonlinearElasticityULMixed.h | 3 ++- Apps/FiniteDefElasticity/SIMFiniteDefEl.C | 2 +- Apps/FiniteDefElasticity/accKM.f | 15 +++++++++--- Apps/FiniteDefElasticity/accKMx.f | 2 ++ 5 files changed, 36 insertions(+), 10 deletions(-) diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C index 1044474f..df2ef46a 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C @@ -29,8 +29,8 @@ extern "C" { const double* Sig, double* D, const int& ipsw, const int& iwr); //! \brief Accumulates material stiffness contributions for 2D problems. - void acckm2d_(const int& nEN, const double* dNdx, - const double* D, double* eKt); + void acckm2d_(const int& axS, const int& nEN, const double* Nr, + const double* dNdx, const double* D, double* eKt); //! \brief Accumulates material stiffness contributions for 3D problems. void acckm3d_(const int& nEN, const double* dNdx, const double* D, double* eKt); @@ -193,8 +193,9 @@ const Vector& NonlinearElasticityULMixed::MixedElmMats::getRHSVector () const } -NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n) - : NonlinearElasticityUL(n), Fbar(3), Dmat(7,7) +NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n, + bool axS) + : NonlinearElasticityUL(n,axS), Fbar(3), Dmat(7,7) { if (myMats) delete myMats; myMats = new MixedElmMats(); @@ -422,13 +423,18 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, this->formKG(*eKg,fe.N1,dNdx,r,Sigma,dVol); } + // Radial shape function for axi-symmetric problems + Vector Nr(fe.N1); + if (axiSymmetry && r > 0.0) + Nr /= r; + #ifdef USE_FTNMAT mdma3d_(Bpres,Mpres,Sig.ptr(),Dmat.ptr(),INT_DEBUG,6); // Integrate the material stiffness matrix Dmat *= dVol; if (nsd == 2) - acckm2d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); + acckm2d_(axiSymmetry,Nr.size(),Nr.ptr(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); else acckm3d_(fe.N1.size(),dNdx.ptr(),Dmat.ptr(),eKm->ptr()); #endif @@ -441,6 +447,7 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, Dmat(7,7) /= (Theta*Theta); for (a = 1; a <= fe.N1.size(); a++) for (b = 1; b <= fe.N2.size(); b++) + { for (i = 1; i <= nsd; i++) { myMats->A[Kut](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * Dmat(i,7); @@ -455,7 +462,14 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt, myMats->A[Kut](nsd*(a-1)+k,b) += dNdx(a,4-k)*fe.N2(b) * Dmat(6,7); } myMats->A[Kup](nsd*(a-1)+i,b) += dNdx(a,i)*fe.N2(b) * dVup; + } + if (axiSymmetry) + { + double Nr = fe.N1(a)/r; + myMats->A[Kut](nsd*(a-1)+1,b) += Nr*fe.N2(b) * Dmat(3,7); + myMats->A[Kup](nsd*(a-1)+1,b) += Nr*fe.N2(b) * dVup; } + } myMats->A[Ktt].outer_product(fe.N2,fe.N2*Dmat(7,7),true); // += N2*N2^T*D77 myMats->A[Ktp].outer_product(fe.N2,fe.N2*(-detJW),true); // -= N2*N2^T*|J| diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h index fd193edb..dd6f18af 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.h @@ -45,7 +45,8 @@ class NonlinearElasticityULMixed : public NonlinearElasticityUL public: //! \brief The default constructor invokes the parent class constructor. //! \param[in] n Number of spatial dimensions - NonlinearElasticityULMixed(unsigned short int n = 3); + //! \param[in] axS \e If \e true, and axisymmetric 3D formulation is assumed + NonlinearElasticityULMixed(unsigned short int n = 3, bool axS = false); //! \brief Empty destructor. virtual ~NonlinearElasticityULMixed() {} diff --git a/Apps/FiniteDefElasticity/SIMFiniteDefEl.C b/Apps/FiniteDefElasticity/SIMFiniteDefEl.C index 6e69110a..28a3a42c 100644 --- a/Apps/FiniteDefElasticity/SIMFiniteDefEl.C +++ b/Apps/FiniteDefElasticity/SIMFiniteDefEl.C @@ -38,7 +38,7 @@ SIMFiniteDefEl2D::SIMFiniteDefEl2D (const std::vector& options) break; case SIM::MIXED_QnQn1: nf[1] = 2; // continuous volumetric change and pressure fields - myProblem = new NonlinearElasticityULMixed(2); + myProblem = new NonlinearElasticityULMixed(2,axiSymmetry); break; case SIM::MIXED_QnPn1: // Local discontinuous volumetric change and pressure fields diff --git a/Apps/FiniteDefElasticity/accKM.f b/Apps/FiniteDefElasticity/accKM.f index 9a5c7325..a0eb2f5a 100644 --- a/Apps/FiniteDefElasticity/accKM.f +++ b/Apps/FiniteDefElasticity/accKM.f @@ -11,12 +11,12 @@ C! \brief Material stiffness accumulation routines. C! C======================================================================= C - subroutine accKM2D (NENOD,Shpf,D,EKt) + subroutine accKM2D (axiSymm,NENOD,Shpr,Shpf,D,EKt) C implicit none C - integer NENOD - double precision Shpf(NENOD,2), D(7,7) + integer axiSymm, NENOD + double precision Shpr(NENOD), Shpf(NENOD,2), D(7,7) double precision EKt(2,NENOD,2,NENOD) C integer a, b, i, j @@ -30,6 +30,10 @@ C BBC(2,j) = Shpf(a,1)*D(4,j) & + Shpf(a,2)*D(2,j) end do +C + if (axiSymm .eq. 1) then + BBC(1,:) = BBC(1,:) + Shpr(a)*D(3,1:4) + end if C do b = 1, NENOD C @@ -48,6 +52,11 @@ C EKt(2,a,2,b) = EKt(2,a,2,b) & + BBC(2,4) * Shpf(b,1) & + BBC(2,2) * Shpf(b,2) +C + if (axiSymm .eq. 1) then + EKt(1,a,1,b) = EKt(1,a,1,b) + BBC(1,3)*Shpr(b) + EKt(2,a,1,b) = EKt(2,a,1,b) + BBC(2,3)*Shpr(b) + end if C end do C diff --git a/Apps/FiniteDefElasticity/accKMx.f b/Apps/FiniteDefElasticity/accKMx.f index ce1c21cb..c7265258 100644 --- a/Apps/FiniteDefElasticity/accKMx.f +++ b/Apps/FiniteDefElasticity/accKMx.f @@ -33,6 +33,7 @@ C & + Shpf(a,2)*D(2,j) & + Shpbar(a,2)*D(7,j) end do +C if (axiSymm .eq. 1) then BBC(1,:) = BBC(1,:) + Shpr(a)*D(3,:) end if @@ -63,6 +64,7 @@ C EKt(1,a,1,b) = EKt(1,a,1,b) + BBC(1,3)*Shpr(b) EKt(2,a,1,b) = EKt(2,a,1,b) + BBC(2,3)*Shpr(b) end if +C end do C end do