esys.modellib.mechanics Package¶

Classes¶

DruckerPrager
IterationDivergenceError
LinearPDE
Mechanics
Model

class esys.modellib.mechanics.DruckerPrager(**kwargs)¶

__init__(**kwargs)¶: set up model

doInitialization()¶

doStep(dt)¶

doStepPostprocessing(dt)¶: accept all the values:

doStepPreprocessing(dt)¶

setStress()¶

setTangentialTensor()¶

class esys.modellib.mechanics.IterationDivergenceError¶

Exception which is thrown if there is no convergence of the iteration process at a time step.

But there is a chance that a smaller step could help to reach convergence.

__init__(*args, **kwargs)¶

class esys.modellib.mechanics.LinearPDE(domain, numEquations=None, numSolutions=None, isComplex=False, debug=False)¶

This class is used to define a general linear, steady, second order PDE for an unknown function u on a given domain defined through a Domain object.

For a single PDE having a solution with a single component the linear PDE is defined in the following form:

-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)

where grad(F) denotes the spatial derivative of F. Einstein’s summation convention, ie. summation over indexes appearing twice in a term of a sum performed, is used. The coefficients A, B, C, D, X and Y have to be specified through Data objects in Function and the coefficients A_reduced, B_reduced, C_reduced, D_reduced, X_reduced and Y_reduced have to be specified through Data objects in ReducedFunction. It is also allowed to use objects that can be converted into such Data objects. A and A_reduced are rank two, B, C, X, B_reduced, C_reduced and X_reduced are rank one and D, D_reduced, Y and Y_reduced are scalar.

The following natural boundary conditions are considered:

n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y

where n is the outer normal field. Notice that the coefficients A, A_reduced, B, B_reduced, X and X_reduced are defined in the PDE. The coefficients d and y are each a scalar in FunctionOnBoundary and the coefficients d_reduced and y_reduced are each a scalar in ReducedFunctionOnBoundary.

Constraints for the solution prescribe the value of the solution at certain locations in the domain. They have the form

u=r where q>0

r and q are each scalar where q is the characteristic function defining where the constraint is applied. The constraints override any other condition set by the PDE or the boundary condition.

The PDE is symmetrical if

A[i,j]=A[j,i] and B[j]=C[j] and A_reduced[i,j]=A_reduced[j,i] and B_reduced[j]=C_reduced[j]

For a system of PDEs and a solution with several components the PDE has the form

-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i]

A and A_reduced are of rank four, B, B_reduced, C and C_reduced are each of rank three, D, D_reduced, X_reduced and X are each of rank two and Y and Y_reduced are of rank one. The natural boundary conditions take the form:

n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]

The coefficient d is of rank two and y is of rank one both in FunctionOnBoundary. The coefficients d_reduced is of rank two and y_reduced is of rank one both in ReducedFunctionOnBoundary.

Constraints take the form

u[i]=r[i] where q[i]>0

r and q are each rank one. Notice that at some locations not necessarily all components must have a constraint.

The system of PDEs is symmetrical if

A[i,j,k,l]=A[k,l,i,j]

A_reduced[i,j,k,l]=A_reduced[k,l,i,j]

B[i,j,k]=C[k,i,j]

B_reduced[i,j,k]=C_reduced[k,i,j]

D[i,k]=D[i,k]

D_reduced[i,k]=D_reduced[i,k]

d[i,k]=d[k,i]

d_reduced[i,k]=d_reduced[k,i]

LinearPDE also supports solution discontinuities over a contact region in the domain. To specify the conditions across the discontinuity we are using the generalised flux J which, in the case of a system of PDEs and several components of the solution, is defined as

J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]

For the case of single solution component and single PDE J is defined as

J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]

In the context of discontinuities n denotes the normal on the discontinuity pointing from side 0 towards side 1 calculated from FunctionSpace.getNormal of FunctionOnContactZero. For a system of PDEs the contact condition takes the form

n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]

where J0 and J1 are the fluxes on side 0 and side 1 of the discontinuity, respectively. jump(u), which is the difference of the solution at side 1 and at side 0, denotes the jump of u across discontinuity along the normal calculated by jump. The coefficient d_contact is of rank two and y_contact is of rank one both in FunctionOnContactZero or FunctionOnContactOne. The coefficient d_contact_reduced is of rank two and y_contact_reduced is of rank one both in ReducedFunctionOnContactZero or ReducedFunctionOnContactOne. In case of a single PDE and a single component solution the contact condition takes the form

n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)

In this case the coefficient d_contact and y_contact are each scalar both in FunctionOnContactZero or FunctionOnContactOne and the coefficient d_contact_reduced and y_contact_reduced are each scalar both in ReducedFunctionOnContactZero or ReducedFunctionOnContactOne.

Typical usage:

p = LinearPDE(dom)
p.setValue(A=kronecker(dom), D=1, Y=0.5)
u = p.getSolution()

__init__(domain, numEquations=None, numSolutions=None, isComplex=False, debug=False)¶

Initializes a new linear PDE.

Parameters

domain (Domain) – domain of the PDE
numEquations – number of equations. If None the number of equations is extracted from the PDE coefficients.
numSolutions – number of solution components. If None the number of solution components is extracted from the PDE coefficients.
debug – if True debug information is printed

checkSymmetry(verbose=True)¶

Tests the PDE for symmetry.

Parameters: verbose (bool) – if set to True or not present a report on coefficients which break the symmetry is printed.
Returns: True if the PDE is symmetric
Return type: bool
Note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.

createOperator()¶: Returns an instance of a new operator.

getFlux(u=None)¶

Returns the flux J for a given u.

J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]

or

J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]

Parameters: u (Data or None) – argument in the flux. If u is not present or equals None the current solution is used.
Returns: flux
Return type: Data

getRequiredOperatorType()¶: Returns the system type which needs to be used by the current set up.

getResidual(u=None)¶

Returns the residual of u or the current solution if u is not present.

Parameters: u (Data or None) – argument in the residual calculation. It must be representable in self.getFunctionSpaceForSolution(). If u is not present or equals None the current solution is used.
Returns: residual of u
Return type: Data

getSolution()¶

Returns the solution of the PDE.

Returns: the solution
Return type: Data

getSystem()¶

Returns the operator and right hand side of the PDE.

Returns: the discrete version of the PDE
Return type: tuple of Operator and Data

insertConstraint(rhs_only=False)¶

Applies the constraints defined by q and r to the PDE.

Parameters: rhs_only (bool) – if True only the right hand side is altered by the constraint

setValue(**coefficients)¶

Sets new values to coefficients.

Parameters

coefficients – new values assigned to coefficients
A (any type that can be cast to a Data object on Function) – value for coefficient A
A_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient A_reduced
B (any type that can be cast to a Data object on Function) – value for coefficient B
B_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient B_reduced
C (any type that can be cast to a Data object on Function) – value for coefficient C
C_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient C_reduced
D (any type that can be cast to a Data object on Function) – value for coefficient D
D_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient D_reduced
X (any type that can be cast to a Data object on Function) – value for coefficient X
X_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient X_reduced
Y (any type that can be cast to a Data object on Function) – value for coefficient Y
Y_reduced (any type that can be cast to a Data object on ReducedFunction) – value for coefficient Y_reduced
d (any type that can be cast to a Data object on FunctionOnBoundary) – value for coefficient d
d_reduced (any type that can be cast to a Data object on ReducedFunctionOnBoundary) – value for coefficient d_reduced
y (any type that can be cast to a Data object on FunctionOnBoundary) – value for coefficient y
d_contact (any type that can be cast to a Data object on FunctionOnContactOne or FunctionOnContactZero) – value for coefficient d_contact
d_contact_reduced (any type that can be cast to a Data object on ReducedFunctionOnContactOne or ReducedFunctionOnContactZero) – value for coefficient d_contact_reduced
y_contact (any type that can be cast to a Data object on FunctionOnContactOne or FunctionOnContactZero) – value for coefficient y_contact
y_contact_reduced (any type that can be cast to a Data object on ReducedFunctionOnContactOne or ReducedFunctionOnContactZero) – value for coefficient y_contact_reduced
d_dirac (any type that can be cast to a Data object on DiracDeltaFunctions) – value for coefficient d_dirac
y_dirac (any type that can be cast to a Data object on DiracDeltaFunctions) – value for coefficient y_dirac
r (any type that can be cast to a Data object on Solution or ReducedSolution depending on whether reduced order is used for the solution) – values prescribed to the solution at the locations of constraints
q (any type that can be cast to a Data object on Solution or ReducedSolution depending on whether reduced order is used for the representation of the equation) – mask for location of constraints

Raises

IllegalCoefficient – if an unknown coefficient keyword is used

class esys.modellib.mechanics.Mechanics(**kwargs)¶

base class for mechanics models in updated lagrangean framework

Note: Instance variable domain - domain (in)
Note: Instance variable internal_force - =Data()
Note: Instance variable external_force - =Data()
Note: Instance variable prescribed_velocity - =Data()
Note: Instance variable location_prescribed_velocity - =Data()
Note: Instance variable temperature - = None
Note: Instance variable expansion_coefficient - = 0.
Note: Instance variable bulk_modulus - =1.
Note: Instance variable shear_modulus - =1.
Note: Instance variable rel_tol - =1.e-3
Note: Instance variable abs_tol - =1.e-15
Note: Instance variable max_iter - =10
Note: Instance variable displacement - =None
Note: Instance variable stress - =None

__init__(**kwargs)¶

set up the model

Parameters: debug (bool) – debug flag

SAFTY_FACTOR_ITERATION = 0.01¶

doInitialization()¶: initialize model

doStep(dt)¶

doStepPostprocessing(dt)¶: accept all the values:

doStepPreprocessing(dt)¶

step up pressure iteration

if run within a time dependend problem extrapolation of pressure from previous time steps is used to get an initial guess (that needs some work!!!!!!!)

getSafeTimeStepSize(dt)¶: returns new step size

terminateIteration()¶: iteration is terminateIterationd if relative pressure change is less than rel_tol

class esys.modellib.mechanics.Model(parameters=[], **kwargs)¶

A Model object represents a process marching over time until a finalizing condition is fulfilled. At each time step an iterative process can be performed and the time step size can be controlled. A Model has the following work flow:

doInitialization()
while not terminateInitialIteration(): doInitialStep()
doInitialPostprocessing()
while not finalize():
    dt=getSafeTimeStepSize(dt)
    doStepPreprocessing(dt)
    while not terminateIteration(): doStep(dt)
    doStepPostprocessing(dt)
doFinalization()

where doInitialization, finalize, getSafeTimeStepSize, doStepPreprocessing, terminateIteration, doStepPostprocessing, doFinalization are methods of the particular instance of a Model. The default implementations of these methods have to be overwritten by the subclass implementing a Model.

__init__(parameters=[], **kwargs)¶

Creates a model.

Just calls the parent constructor.

UNDEF_DT = 1e+300¶

doFinalization()¶

Finalizes the time stepping.