..and finally the Qq/Qq-1 formulations

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1273 e10b68d5-8a6e-419e-a041-bce267b0401d

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,
+  void acckm2d_(const int& axS, const int& nEN, const double* Nr,
-		const double* D, double* eKt);
+		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)
+NonlinearElasticityULMixed::NonlinearElasticityULMixed (unsigned short int n,
-  : NonlinearElasticityUL(n), Fbar(3), Dmat(7,7)
+							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);
@@ -456,6 +463,13 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
 	}
 	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|

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() {}
 

Apps/FiniteDefElasticity/SIMFiniteDefEl.C

@@ -38,7 +38,7 @@ SIMFiniteDefEl2D::SIMFiniteDefEl2D (const std::vector<int>& 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

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
+      integer          axiSymm, NENOD
-      double precision Shpf(NENOD,2), D(7,7)
+      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

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