[CIG-SHORT] Plane strain test problem

[Brad Aagaard]

Birendra,

The main problem is that you have Dirichlet BC for z only. Without any 
constraints in the x and y directions, the solution can rotate or 
translate in the xy plane with rigid body motion. Rather than two 
Neumann BC that apply equal and opposite tractions, replace one with a 
Dirichlet BC to remove rigid body motion in the y direction and add a an 
additional Dirichlet BC to remove rigid body motion in the x direction 
(may only need to pin a couple nodes).

A minor issue is that your z_neg and z_pos node sets do not cover the 
entire face. This means that not all z DOF are fixed so it is not a true 
plane strain problem.

Brad


On 05/11/2012 03:58 PM, Birendra jha wrote:
> Dear Brad,
>
> I am running this plane strain test case (files attached) with equal
> compression on two y faces, zero tractions on x faces, and zero
> normal displacement on z faces. (Does this replicate Plane Strain?)
>
> I don't understand why the y-displacements are not uniform along the
> length i.e. along x. Why does the grid rotate counterclockwise when
> the compressions are equal on opposite sides?
>
> Any help will be greatly appreciated.
>
> Thanks and best regards Birendra
>
> --- On Thu, 5/10/12, Brad Aagaard<baagaard at usgs.gov>  wrote:
>
>> From: Brad Aagaard<baagaard at usgs.gov> Subject: Re: Linear system
>> solve To: "Birendra jha"<bjha7333 at yahoo.com> Cc:
>> cig-short at geodynamics.org Date: Thursday, May 10, 2012, 2:22 AM On
>> 05/09/2012 01:23 PM, Birendra jha wrote:
>>> Dear developers,
>>>
>>> I have a question about KSPSolve in,
>>>
>>> err = KSPSolve(_ksp, residualVec, solutionVec);
>>> CHECK_PETSC_ERROR(err);
>>>
>>> called inside,
>>>
>>> pylith::problems::SolverLinear::solve(
>>> topology::Field<topology::Mesh>* solution,
>> topology::Jacobian*
>>> jacobian, const
>> topology::Field<topology::Mesh>&   residual)
>>>
>>> I find that for same residualVec and jacobian matrix,
>> my solutionVec
>>> (i.e. dispIncr) changes. What could be the reason? Does
>> KSPSolve do
>>> anything other than solving Ax=b system, after it's
>> provided with A
>>> and b? I checked the jacobian in Implicit.py:step() using
>>> self.jacobian.view(), before the call to solve() is
>> made. I checked
>>> residualVec and solutionVec using
>>> VecView(solutionVec,PETSC_VIEWER_STDOUT_WORLD) inside
>> the
>>> SolverLinear::solve().
>>>
>>
>> What do you mean solutionVec changes? Do you mean the address of
>> solutionVec changes or values inside the solutionVec change? When
>> do they change?
>>
>> KSPSolve is a PETSc routine. The documentation is here
>> http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/KSP/KSPSolve.html
>>
>>
>>
Within SolverLinear::Solve() we do have to transfer values
>> between Sieve sections and PETSc vectors via
>> scatterSectionToVector() and scatterVectorToSection(). These calls
>> will change values of the corresponding Sieve sections or PETSc
>> vectors.
>>
>>> Boundary conditions are independent of time.
>> Compression on two side
>>> boundaries, and nothing specified for other two sides
>> which means
>>> traction-free surface, is that correct?
>>
>> If not BC is specified on a surface, it will be traction free as in
>> nearly all finite-element implementations.
>>
>> Brad
>>
>>