]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
step:44 fixed intial J projection error, updated documentation, intro and results
authormcbride <mcbride@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 24 Feb 2012 16:26:00 +0000 (16:26 +0000)
committermcbride <mcbride@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 24 Feb 2012 16:26:00 +0000 (16:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@25164 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-44/doc/intro.dox
deal.II/examples/step-44/doc/results.dox
deal.II/examples/step-44/step-44.cc

index 26bfd351110a32103a1e6cc92d5d6f28018e608b..6c4a51a43d6ca5722438ef70c03bcabfa1b0bbba 100644 (file)
@@ -274,7 +274,8 @@ The Helmholtz free energy corresponding to a compressible <a href="http://en.wik
 @f]
 where $\kappa := \lambda + 2/3 \mu$ is the bulk modulus ($\lambda$ and $\mu$ are the Lame parameters)
 and $\overline{I}_1 := \textrm{tr}\ \overline{\mathbf{b}}$.
-The function $\mathcal{G}(J)$ is required to be strictly convex and satisfy the condition $\mathcal{G}(1) = 0$.
+The function $\mathcal{G}(J)$ is required to be strictly convex and satisfy the condition $\mathcal{G}(1) = 0$, 
+among others, see Holzapfel (2001) for further details.
 In this work $\mathcal{G}:=\frac{1}{4} [ J^2 - 1 - 2\textrm{ln}J ]$.
 
 Incompressibility imposes the isochoric constraint that $J=1$ for all motions $\mathbf{\varphi}$.
@@ -430,8 +431,14 @@ The third is the definition of the pressure $\widetilde{p}$.
 &=
 - \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V 
 + \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V
- + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \,.
+ + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \\
+&=
+- \int_{\Omega_0} [\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} :\boldsymbol{\tau}~\mathrm{d}V 
++ \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V
+ + \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \, ,
 @f}
+where 
+$[\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} = 1/2[ \mathrm{grad}\delta\mathbf{u} + [\mathrm{grad}\delta\mathbf{u}]^T] $.
 
 We will use an iterative Newton-Raphson method to solve the nonlinear residual equation $R$.
 For the sake of simplicity we assume dead loading, i.e. the loading does not change due to the deformation.
index 4198d3195062b58e8b21650a5e3f95b5ccf87163..7726047dec1b522e8ce8be127f3d4dff01bde08b 100644 (file)
@@ -34,7 +34,7 @@ as it should
 A typical screen output generated by running the problem is shown below. 
 The particular case demonstrated is that of the $Q_2-DGPM_1-DGPM_1$ formulation.
 It is clear that, using the Newton-Raphson method, quadratic convergence of the solution is obtained. 
-Other than in the first time step, solution convergence is achieved within 5 Newton increments. 
+Solution convergence is achieved within 5 Newton increments for all time-steps. 
 The converged displacement's $L_2$-norm is several orders of magnitude less than the geometry scale.
 
 @code
@@ -49,18 +49,15 @@ Timestep 1 @ 0.1s
 ___________________________________________________________________________________________________________________________________________________________
                  SOLVER STEP                   |  LIN_IT   LIN_RES    RES_NORM     RES_U     RES_P      RES_J     NU_NORM      NU_U       NU_P       NU_J 
 ___________________________________________________________________________________________________________________________________________________________
-  0  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     516  2.432e-06  1.000e+00  1.000e+00  0.000e+00  0.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  
-  1  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     290  4.660e-03  5.446e+00  8.554e-02  2.602e-10  1.437e+01  4.191e+00  2.550e+00  4.191e+00  1.456e+04  
-  2  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     339  7.754e-03  2.780e+00  1.901e+00  7.310e-12  5.352e+00  3.597e+00  2.129e+00  3.597e+00  1.316e+03  
-  3  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     325  2.931e-03  1.443e+00  1.380e+00  4.568e-12  1.115e+00  9.874e-01  4.266e-01  9.874e-01  1.250e+03  
-  4  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     390  2.072e-04  8.894e-02  8.523e-02  2.469e-13  6.711e-02  4.998e-02  2.656e-02  4.998e-02  1.237e+02  
-  5  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     506  8.187e-07  4.081e-04  3.921e-04  8.575e-16  2.980e-04  2.522e-04  1.344e-04  2.522e-04  5.493e-01  
-  6  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     613  1.665e-11  3.105e-07  3.105e-07  2.613e-20  5.834e-09  3.108e-07  3.428e-08  3.108e-07  1.066e-05  
-  7  ASM_R  CONVERGED! 
+  0  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     786  2.118e-06  1.000e+00  1.000e+00  0.000e+00  0.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  
+  1  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     552  1.031e-03  8.563e-02  8.563e-02  9.200e-13  3.929e-08  1.060e-01  3.816e-02  1.060e-01  1.060e-01  
+  2  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     667  5.602e-06  2.482e-03  2.482e-03  3.373e-15  2.982e-10  2.936e-03  2.053e-04  2.936e-03  2.936e-03  
+  3  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     856  6.469e-10  2.129e-06  2.129e-06  2.245e-19  1.244e-13  1.887e-06  7.289e-07  1.887e-06  1.887e-06  
+  4  ASM_R  CONVERGED! 
 ___________________________________________________________________________________________________________________________________________________________
 Relative errors:
-Displacement:  3.428e-08
-Force:                 6.310e-12
+Displacement:  7.289e-07
+Force:                 2.451e-10
 Dilatation:    1.353e-07
 v / V_0:       1.000e-09 / 1.000e-09 = 1.000e+00
 
@@ -71,17 +68,17 @@ Timestep 10 @ 1.000e+00s
 ___________________________________________________________________________________________________________________________________________________________
                  SOLVER STEP                   |  LIN_IT   LIN_RES    RES_NORM     RES_U     RES_P      RES_J     NU_NORM      NU_U       NU_P       NU_J 
 ___________________________________________________________________________________________________________________________________________________________
-  0  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     556  2.542e-06  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  
-  1  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     455  8.017e-05  1.538e-01  1.538e-01  3.522e+11  4.617e+07  6.010e-02  7.396e-02  6.010e-02  6.011e-02  
-  2  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     554  4.012e-07  2.896e-03  2.896e-03  2.127e+09  5.117e+04  2.720e-03  1.430e-03  2.720e-03  2.724e-03  
-  3  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     620  4.946e-10  1.728e-06  1.728e-06  1.315e+06  9.778e+01  1.150e-06  6.268e-07  1.150e-06  1.150e-06  
+  0  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     874  2.358e-06  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  1.000e+00  
+  1  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     658  2.942e-04  1.544e-01  1.544e-01  1.208e+13  1.855e+06  6.014e-02  7.398e-02  6.014e-02  6.014e-02  
+  2  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     790  2.206e-06  2.908e-03  2.908e-03  7.302e+10  2.067e+03  2.716e-03  1.433e-03  2.716e-03  2.717e-03  
+  3  ASM_R  ASM_K  CST  ASM_SC  SLV  PP  UQPH  |     893  2.374e-09  1.919e-06  1.919e-06  4.527e+07  4.100e+00  1.672e-06  6.842e-07  1.672e-06  1.672e-06  
   4  ASM_R  CONVERGED! 
 ___________________________________________________________________________________________________________________________________________________________
 Relative errors:
-Displacement:  6.268e-07
-Force:                 1.873e-10
-Dilatation:    1.525e-06
-v / V_0:       9.999e-10 / 1.000e-09 = 9.999e-01
+Displacement:  6.842e-07
+Force:                 8.995e-10
+Dilatation:    1.528e-06
+v / V_0:       1.000e-09 / 1.000e-09 = 1.000e+00
 @endcode
 
 
@@ -229,7 +226,7 @@ However, it is the author's experience that the Jacobi method of preconditioning
      <td align="center">
         @image html "step-44.Normalised_runtime.png"
        <p align="center">
-        Runtime on a 4-core machine, normalised against the lowest grid resolution Q1-P0 solution that utilised a SSOR preconditioner.
+        Runtime on a 4-core machine, normalised against the lowest grid resolution $Q_1-DGPM_0-DGPM_0$ solution that utilised a SSOR preconditioner.
        </p>
     </td>
   </tr>
index b61ff2dda373a4b79840d49ef3617077d174be7c..6745bf7075ed30b932eceee8bbe11bd242d179f4 100644 (file)
@@ -1363,14 +1363,45 @@ namespace Step44
 // interchangeable. We choose to increment time linearly using a constant time
 // step size.
 //
-// We start the function with preprocessing, and then output the initial grid
-// before starting the simulation proper with the first time (and loading)
-// increment:
+// We start the function with preprocessing, setting the initial dilatation
+// values, and then output the initial grid before starting the simulation
+//  proper with the first time (and loading)
+// increment.
+//
+// Care must be taken (or at least some thought given) when imposing the
+// constraint $\widetilde{J}=1$ on the initial solution field. The constraint
+// corresponds to the determinant of the deformation gradient in the undeformed
+// configuration, which is the identity tensor.
+// We use FE_DGPMonomial bases to interpolate the dilatation field, thus we can't
+// simply set the corresponding dof to unity as they correspond to the
+// monomial coefficients. Thus we use the VectorTools::project function to do
+// the work for us. The VectorTools::project function requires an argument
+// indicating the hanging node constraints. We have none in this program
+// So we have to create a constraint object. In its original state, constraint
+// objects are unsorted, and have to be sorted (using the ConstraintMatrix::close function)
+// before they can be used. Have a look at step-21 for more information.
+// We only need to enforce the initial condition on the dilatation.
+// In order to do this, we make use of a ComponentSelectFunction which acts
+// as a mask and sets the J_component of n_components to 1. This is exactly what
+// we want. Have a look at its usage in step-20 for more information.
   template <int dim>
   void Solid<dim>::run()
   {
     make_grid();
     system_setup();
+    {
+      ConstraintMatrix constraints;
+      constraints.close();
+
+      const ComponentSelectFunction<dim>
+        J_mask (J_component, n_components);
+
+      VectorTools::project (dof_handler_ref,
+                            constraints,
+                            QGauss<dim>(degree+2),
+                            J_mask,
+                            solution_n);
+    }
     output_results();
     time.increment();
 
@@ -1647,6 +1678,9 @@ namespace Step44
 // this operation (we could, in principle simply create a new task using
 // Threads::new_task for each cell) but there is not much harm done to doing
 // it this way anyway.
+// Furthermore, should their be different material models associated with a
+// quadrature point, requiring varying levels of computational expense, then
+// the method used here could be advantageous.
   template <int dim>
   struct Solid<dim>::PerTaskData_UQPH
   {
@@ -1851,16 +1885,12 @@ namespace Step44
 
     tangent_matrix.reinit(sparsity_pattern);
 
-                                    // We then set up storage vectors noting
-                                    // that the dilatation is unity
-                                    // (i.e. $\widetilde{J} = 1$) in the
-                                    // undeformed configuration...
+                                    // We then set up storage vectors
     system_rhs.reinit(dofs_per_block);
     system_rhs.collect_sizes();
 
     solution_n.reinit(dofs_per_block);
     solution_n.collect_sizes();
-    solution_n.block(J_dof) = 1.0;
 
                                     // ...and finally set up the quadrature
                                     // point history:

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.