Projection of function to mixed variational form error - dolfinx

variational formulation
1

Hi, I’m trying to convert my FEniCS code to FEniCSx. In particular I’m trying to implement a mixed formulation (strain/displacement) using stabilization terms in order to use the same interpolation in both fields (P1/P1). This wasn’t a problem in FEniCS using project() in the variational form, but in FEniCSx I’m getting the following error:

Info    : Reading 'mesh/cooks_membrane_2.msh'...
Info    : 9 nodes
Info    : 12 elements
Info    : Done reading 'mesh/cooks_membrane_2.msh'
Traceback (most recent call last):
  File "/home/batipc/fenicsx/1-TEST_linear/sample.py", line 121, in <module>
    FF = -tau_eps(W)*(as_tensor(gamma[i,j]*C0[i,j,k,l]*Orth_proj(W_eps,sym(nabla_grad(u)),sym(nabla_grad(uk)))[k,l]))*dx \
  File "/home/batipc/fenicsx/1-TEST_linear/sample.py", line 84, in Orth_proj
    proj = (res - proj2space(res_k, W))
  File "/home/batipc/fenicsx/1-TEST_linear/sample.py", line 46, in proj2space
    L  = form(inner(v, w) * dx)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 176, in form
    return _create_form(form)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 171, in _create_form
    return _form(form)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 145, in _form
    ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/jit.py", line 56, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 186, in compile_forms
    impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 250, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/compiler.py", line 102, in compile_ufl_objects
    ir = compute_ir(analysis, object_names, prefix, options, visualise)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/representation.py", line 200, in compute_ir
    irs = [
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/representation.py", line 201, in <listcomp>
    _compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/representation.py", line 485, in _compute_integral_ir
    integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/integral.py", line 72, in compute_integral_ir
    S = build_scalar_graph(expression)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/analysis/graph.py", line 80, in build_scalar_graph
    scalar_expressions = rebuild_with_scalar_subexpressions(G)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/analysis/graph.py", line 125, in rebuild_with_scalar_subexpressions
    V_symbols = value_numberer.compute_symbols()
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/analysis/valuenumbering.py", line 68, in compute_symbols
    symbol = f(expr)
  File "/home/batipc/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/ffcx/ir/analysis/valuenumbering.py", line 184, in _modified_terminal
    raise RuntimeError("Internal error in value numbering.")
RuntimeError: Internal error in value numbering.

In this case, I’m using a projection similar to the one used in dolfiny here https://github.com/michalhabera/dolfiny/blob/master/dolfiny/projection.py. Here’s my minimal example and geometry for the msh file (here I’m using only the stabilization term to reproduce the error):

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
3
1 1 "izq"
1 2 "der"
2 1 "surface"
$EndPhysicalNames
$Nodes
9
1 0 0 0
2 0 44 0
3 48 60 0
4 48 44 0
5 0 22 0
6 23.99999999999994 51.99999999999998 0
7 48 52 0
8 23.99999999999996 21.99999999999996 0
9 23.99999999999995 36.99999999999997 0
$EndNodes
$Elements
12
1 1 2 1 1 1 5
2 1 2 1 1 5 2
3 1 2 2 3 3 7
4 1 2 2 3 7 4
5 2 2 1 6 1 5 8
6 2 2 1 6 8 5 9
7 2 2 1 6 8 9 4
8 2 2 1 6 4 9 7
9 2 2 1 6 5 2 9
10 2 2 1 6 9 2 6
11 2 2 1 6 9 6 7
12 2 2 1 6 7 6 3
$EndElements


import os
import numpy as np
from ufl import (FiniteElement,VectorElement, TensorElement, MixedElement, Measure, TestFunction, TrialFunction,
                 split, as_tensor, dot, inner, nabla_grad,  sym, indices, derivative, CellDiameter, dx, ds)
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.fem import (Constant, Function, FunctionSpace, VectorFunctionSpace, TensorFunctionSpace,
                         assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc,assemble,Expression)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, 
                               create_vector, create_matrix, set_bc)
from dolfinx.io import (XDMFFile, gmshio)

## ***************** Functions *****************
def build_space(_name,_comm,_gdim): 
    #Import mesh
    mesh, _, facets= gmshio.read_from_msh(_name, _comm, gdim=_gdim)
    
    # Build mixed function space  
    P1 = VectorElement("CG",mesh.ufl_cell(),1)                 # displacement
    P2 = TensorElement("CG",mesh.ufl_cell(),1,symmetry = True) # strain

    MX = MixedElement([P1,P2])
    W = FunctionSpace(mesh,MX)

    # Dirichlet boundary conditions
    fdim = mesh.topology.dim - 1
    dofs_i0 = locate_dofs_topological(W.sub(0).sub(0), fdim, facets.find(1))
    dofs_i1 = locate_dofs_topological(W.sub(0).sub(1), fdim, facets.find(1))
    bcs = []
    bc_izq1 = dirichletbc(PETSc.ScalarType((0.0)), dofs_i0, W.sub(0).sub(0)) # clamped side (left) -> Ux(0)=0
    bc_izq2 = dirichletbc(PETSc.ScalarType((0.0)), dofs_i1, W.sub(0).sub(1)) # clamped side (left) -> Uy(0)=0
    bcs.append(bc_izq1);bcs.append(bc_izq2)
    
    return W,bcs,mesh,facets

def proj2space(v, V, bcs=[]):    
    metad = {"quadrature_degree": 2}    
    dx = Measure('dx', domain=mesh, metadata=metad)

    # Define variational problem for projection
    w  = TestFunction(V)
    Pv = TrialFunction(V)
    a  = form(inner(Pv, w) * dx)
    L  = form(inner(v, w) * dx)

    # Assemble linear system
    A = assemble_matrix(a, bcs)
    A.assemble()
    b = assemble_vector(L)
    apply_lifting(b, [a], [bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, bcs)

    # Solve linear system
    solver = PETSc.KSP().create(A.getComm())
    solver.setOperators(A)
    uh  = Function(V)
    solver.solve(b, u.vector)
    return uh

#%% Define constitutive tensor (list)
def VoigtTo4Tensor_2d(A):
        A11, A12, A13 = A[0,0], A[0,1], A[0,2]
        A22, A23  = A[1,1], A[1,2]
        A33 = A[2,2]
        A21 , A31 = A12 , A13
        A32 = A23 
        return  [ [ \
                [ [ A11 , A13 ] , [ A13 , A12 ] ] , \
                [ [ A31 , A33 ] , [ A33 , A32 ] ] , \
                ] , [ \
                [ [ A31 , A33 ] , [ A33 , A32 ] ] , \
                [ [ A21 , A23 ] , [ A23 , A22 ] ] ] ] 

#%% Stabilization parameters
def tau_eps(W):
    h = CellDiameter(W.mesh)
    return ce*h/L0

#%% Orthogonal projection VMS
def Orth_proj(W,res,res_k):
    proj = (res - proj2space(res_k, W))
    return proj

#%% ***************** Main *************************

#import mesh
comm = MPI.COMM_WORLD
name = "mesh/cooks_membrane_2.msh"
gdim = 2
W,bcs,mesh,facets = build_space(name,comm,gdim)  # Facets: 1->left, 2->right

# Scalar variables
E  = Constant(mesh,PETSc.ScalarType(200e3))  # Young's Moduli [MPa]
nu = Constant(mesh,PETSc.ScalarType(0.3))    # Poisson Coefficient
L0 = Constant(mesh,PETSc.ScalarType(50.0))   # Characteristic lenght [mm]
ce = Constant(mesh,PETSc.ScalarType(1.0))    # Arbitrary constant for tau_eps

# Compute constitutive tensor for plain strain
C_voigt = E/(1+nu)/(1-2*nu)*np.array([ [1-nu,  nu , 0],\
                                        [ nu , 1-nu, 0],\
                                        [  0 ,  0  , (1-2*nu)/2]])
C0 = as_tensor(VoigtTo4Tensor_2d(C_voigt))  # Constitutive tensor

#Functions over the mesh
w = Function(W)             # mixed function (u,eps)
u,eps = split(w)            # split mixed function
z = TestFunction(W)         # test function (v,gamma)
v,gamma = split(z)          # split test function
wk = Function(W)            # mixed function (uk,epsk) at previous iteration
uk,epsk = split(wk)         # split mixed function at previous iteration

# Defining variational problem
W_eps,_ = W.sub(1).collapse()
T = Constant(mesh, PETSc.ScalarType((0, 1000))) # Traction forces [N/mm]
i,j,k,l = indices(4)

# First error
ds = Measure('ds', domain=mesh, subdomain_data=facets)
FF = -tau_eps(W)*(as_tensor(gamma[i,j]*C0[i,j,k,l]*Orth_proj(W_eps,sym(nabla_grad(u)),sym(nabla_grad(uk)))[k,l]))*dx \
        - dot(v,T)*ds(2)    # Variational form 

# Second error
# FF2 = -tau_eps(W)*(as_tensor(gamma[i,j]*C0[i,j,k,l]*sym(nabla_grad(u))[k,l]))*dx - dot(v,T)*ds(2)    # Variational form

# JJ = derivative(FF2, w)      # Jacobian
# residual = form(FF2)
# jacob = form(JJ)

# ... newton solver

I’m getting the same error when using instead a mixed variational form that doesn’t use a projection (FF2). Seems like the error comes when using the function form() with the mixed form. Is it not compatible? I’m interested in using this to use a custom newton solver like the one showed in this tutorial: Custom Newton solvers — FEniCSx tutorial. I’m using Dolfinx v0.6.0 in conda

1 Like
2

I’ve made an issue [BUG]: Symmetric spaces broken · Issue #2750 · FEniCS/dolfinx · GitHub with a simplified problem.

1 Like
3

Hello @dokken, and thank you for your repply. I’ve been following the state of the bug and I saw there’s people working on it, so I’m looking forward for the new version release. By the moment I’ve been working on FEniCS Legacy in the implementation of the mixed form with a explicit damage constituitive law, but I ran into a problem when computing this variable at nodal level. The usual procedure in FEM is to compute the constitutive relation at integration points, but the idea of this implementation is to take advantage of the P1P1 formulation so the strain field is continuos, so I should be able to calculate the damage variable (d) at nodes, but the issue comes when using the variable d as a function in a P1 FunctionSpace. When I do that, I get wrong values (d should be between 0 and 1) and the only solution that has worked is to use d in a DP0 FunctionSpace. I think the problem emerge because the projection I’m doing with both internal variables (d and r), but the code doesn’t works if I don’t do those projections. Is there an alternative so I can interpolate instead of project? Here is my example with a two elements geometry (I’m sorry but I didn’t figure out how to make a shorter example):

from dolfin import *
import os
import numpy as np
from ufl import indices,nabla_div

def Meshconv(nombre,geo_file_content):
    #import geometry from GMSH and convertz to xml
    
    try: 
        with open("%s.mesh/%s.geo"%(nombre,nombre), 'w') as filed:
            filed.write(geo_file_content)
    except Exception:
        os.mkdir("%s.mesh"%(nombre))
        with open("%s.mesh/%s.geo"%(nombre,nombre), 'w') as filed:
            filed.write(geo_file_content)   

    os.system('gmsh -2 %s.mesh/%s.geo -format msh2'%(nombre,nombre))
    os.system('dolfin-convert -i gmsh %s.mesh/%s.msh %s.mesh/%s.xml'%(nombre,nombre,nombre,nombre))

    mesh = Mesh("%s.mesh/%s.xml" %(nombre,nombre))
    contorno = MeshFunction("size_t", mesh,"%s.mesh"%(nombre) +"/"+nombre+"_facet_region.xml")
    subdominios= MeshFunction("size_t",mesh,"%s.mesh"%(nombre) +"/"+nombre+"_physical_region.xml")
    
    return mesh,contorno,subdominios

def build_space():
    # Builds the mesh file from gmsh code 
    # Builds the mixed function space and geometrical boundary conditions
    # Returns: W,bcs,mesh,contorno
    # W: Mixed function space
    # bcs: boundary conditions
    # mesh: mesh
    # contorno: facets of the mesh 
    
    geo_file_content = '''
    //Unitary_square
    cl = 1;
    SetFactory("OpenCASCADE");
    s1 = news; Rectangle(s1) = {0,0,0,1,1,0};

    // physical entities
    Physical Curve("borde_izquierdo",1) = {4};
    Physical Curve("borde_derecho",2) = {2};
    Physical Curve("borde_superior",3) = {3};
    Physical Curve("borde_inferior",4) = {1};
    Physical Surface(1) = {s1};

    //Meshing
    Transfinite Curve{2,4} = 2;
    Transfinite Curve{1,3} = 2; 
    Transfinite Surface{s1};

    Mesh 2;

    Mesh.MshFileVersion = 2.2;
    Save StrCat(StrPrefix(General.FileName), ".msh");
    '''
    #Import mesh
    
    nombre = 'Unitary_Square'
    mesh,contorno,subdominio= Meshconv(nombre, geo_file_content)
    
    # Build mixed function space
    
    P1 = VectorElement("P",mesh.ufl_cell(),1) #displacement
    P2 = TensorElement("P",mesh.ufl_cell(),1,symmetry = True) #strain

    MX = MixedElement([P1,P2])
    W = FunctionSpace(mesh,MX)

    # Dirichlet boundary conditions -> related to geometry
    
    bcs = []
    bc_izq = DirichletBC(W.sub(0),Constant((0,0)),contorno,1) #clamped side (left)
    bcs.append(bc_izq)
    
    return W,bcs,mesh,contorno #,subdominio



#%% Initial condition for mixed element
class InitialCondition_w(UserExpression):
    def eval(self, values, x):
        values[0] = 0.0  #ux
        values[1] = 0.0  #uy
        values[2] = 0.0  #ex
        values[3] = 0.0  #exy
        values[4] = 0.0  #eyx
        values[5] = 0.0  #ey
           
    def value_shape(self):
        return (6,)



#%% Define C tensor
def C_tensor(E,nu):
    #Creates a list in the shape of a 4th order tensor, for plain strain
    #then uses as_tensor function to create a tensor object
    
    d = 2
    componentList = []
    for i in range(0,d):
        componentList += [[],]
        for j in range(0,d):
            componentList[i] += [[],]
            for k in range(0,d):
                componentList[i][j] += [[],]
                for l in range(0,d):
                    componentList[i][j][k] += [1,]
    
    # Components of tensor C (plane strain)
    componentList[0][0][0][0] = (E/((1+nu)*(1-2*nu)))*(1-nu)      #xxxx
    componentList[0][0][0][1] = (E/((1+nu)*(1-2*nu)))* 0          #xxxy
    componentList[0][0][1][0] = (E/((1+nu)*(1-2*nu)))* 0          #xxyx
    componentList[0][0][1][1] = (E/((1+nu)*(1-2*nu)))* nu         #xxyy
    componentList[0][1][0][0] = (E/((1+nu)*(1-2*nu)))* 0          #xyxx
    componentList[0][1][0][1] = (E/((1+nu)*(1-2*nu)))* (1-2*nu)/2 #xyxy
    componentList[0][1][1][0] = (E/((1+nu)*(1-2*nu)))* (1-2*nu)/2 #xyyx
    componentList[0][1][1][1] = (E/((1+nu)*(1-2*nu)))* 0          #xyyy
    componentList[1][0][0][0] = (E/((1+nu)*(1-2*nu)))* 0          #yxxx
    componentList[1][0][0][1] = (E/((1+nu)*(1-2*nu)))* (1-2*nu)/2 #yxxy
    componentList[1][0][1][0] = (E/((1+nu)*(1-2*nu)))* (1-2*nu)/2 #yxyx
    componentList[1][0][1][1] = (E/((1+nu)*(1-2*nu)))* 0          #yxyy
    componentList[1][1][0][0] = (E/((1+nu)*(1-2*nu)))* nu         #yyxx
    componentList[1][1][0][1] = (E/((1+nu)*(1-2*nu)))* 0          #yyxy
    componentList[1][1][1][0] = (E/((1+nu)*(1-2*nu)))* 0          #yyyx
    componentList[1][1][1][1] = (E/((1+nu)*(1-2*nu)))*(1-nu)      #yyyy

    return componentList


#%% Define stress
def sigma_3d(e,C):
    # stress for plain strain
    i,j,k,l = indices(4)
    lambda_ = E*nu/(1+nu)/(1-2*nu)
    s = as_tensor(C[i,j,k,l]*e[k,l],(i,j))
    return as_tensor([[s[0,0], s[0,1], 0],
                      [s[0,1], s[1,1], 0],
                      [  0   ,   0   , lambda_*(e[0,0]+e[1,1])]])

def Principal_stress(stress):
    # Get eigenvalues of stress tensor
    # input: projected stress tensor into DP0 or P1 space (if DP0: mesh.num_cells(), if P1: mesh.num_vertices())
    # output: major principal stress (max eigenvalue eig[1] for plain stress, eig[2] for plain strain)
    mesh = stress.function_space().mesh()
    eig = np.linalg.eigvalsh(stress.vector().get_local().reshape([mesh.num_vertices(),3,3]))
    eigen = Function(VectorFunctionSpace(mesh,'P',1,dim=3))
    eigen.vector().set_local(eig.flatten())
    return eigen[2]



#Damage propagation
def my_max(a, b):
    return (a+b+abs(a-b))/2.0

def Macaulay(s1):
    return (s1+abs(s1))/2.

def Rankine(s):
    s1 = Principal_stress(s)
    return Macaulay(s1)

def calc_new_r(eps, r_n,C):
    # Evolution of internal variable r - explicit
    strs = sigma_3d(eps,C)
    stress_eq = Rankine(project(strs,ET1))
    ret = my_max(stress_eq, r_n)
    return ret

def q_dam(r,r0):
    h = CellDiameter(W.mesh())              # Element size
    b = (1-tau_eps(W))*2*h + tau_eps(W)*h   # Bandwidth
    Hd = b/(2*L-b)                          # Softening parameter
    return r0*exp(-2*Hd*((r-r0)/r0))

def calc_d(eps,r_n,C):
    # Damage index evolution 
    r_temp = calc_new_r(eps, r_n,C)
    d_temp = 1.0 - (q_dam(r_temp, r0)/r_temp)
    return d_temp



#%% Stabilization parameters
def tau_eps(W):
    # stabilization parameter strains
    h = CellDiameter(W.mesh())
    return ce*h/L0

def tau_u(W,d):
    # stabilization parameter displacements
    h = CellDiameter(W.mesh())
    return cu*(h*L0)/((1-d)*E)



#%% Operators
def res_e(eps,u):
    # Residual asociated to strain equation
    return (eps-sym(nabla_grad(u)))

def res_u(eps,C):
    # Residual asociated to displacement equation
    i,j,k,l = indices(4)
    return nabla_div(as_tensor(C[i,j,k,l]*eps[k,l],(i,j)))

#%% Orthogonal projection
def Orth_p(W,res,res_k):
    return (res - project(res_k, W))



#%% Define Mixed Variational Problem
def VarForm_full(W,w,z,d,C):  
    # Variational form for incompressible materials
    # Considers all 6 stabilization terms
    (u,eps) = split(w)
    (v,gamma) = split(z)

    T = Constant((0, 0))
    B = Constant((0, 0))

    i,j,k,l = indices(4)

    # ASGS stabilization
    F1 = -as_tensor(gamma[i,j]*C(d)[i,j,k,l]*res_e(eps,u)[k,l])*dx
    stab1 = tau_eps(W)*as_tensor(gamma[i,j]*C(d)[i,j,k,l]*res_e(eps,u)[k,l])*dx
    stab2 = -tau_u(W,d)*inner(nabla_div(as_tensor(C(d)[i,j,k,l]*gamma[k,l],(i,j))),res_u(eps,C(d)))*dx
    F2 = as_tensor(sym(nabla_grad(v))[i,j]*C(d)[i,j,k,l]*eps[k,l])*dx
    stab3 = -tau_eps(W)*as_tensor(sym(nabla_grad(v))[i,j]*C(d)[i,j,k,l]*res_e(eps,u)[k,l])*dx
    F3 = dot(v,B)*dx + dot(v,T)*ds

    FF = F1 + F2 + stab1 + stab2 + stab3 - F3     
    return FF



def solve_damage(W,Escalar,bcs,mesh,contorno,r0,inc_f,file_path):
    
    # Damage index dependent parameters
    C0 = as_tensor(C_tensor(E,nu)) #C tensor assembly (2D)
    C_dam = lambda d : (1-d)*C0   # secant stiffness tensor

    # Functions over the mesh
    w = TrialFunction(W)      # mixed trial function (u,eps)
    w_ = Function(W)          # mixed function (u,eps)
    z = TestFunction(W)       # test function (v,gamma)
    icw_0 = InitialCondition_w(degree=1)     #Initial condition for mixed element
    wk = Function(W)          # mixed function iteration i-1 (u_k,eps_k)
    wk.assign(interpolate(icw_0,W))
    (u_old,eps_old) = wk.split(deepcopy=True) 
    
    ic_0 = Expression('ic',degree=0,ic =0)    #Initial condition
    d = Function(Escalar)                           # d_index (iteration i)
    d.assign(interpolate(ic_0,Escalar))
    dk = Function(Escalar)                          # d_index (iteration i-1)
    dk.assign(interpolate(ic_0,Escalar))

    r = Function(Escalar)                           # r (iteration i)
    r.assign(interpolate(Constant(r0), Escalar))  

    # Global Loop     
    inc = 0
    load.t-=0.001
    while inc <= inc_f:
        
        load.t+= 0.001e-0
        
        L2_error_u = 1.0
        it = 0
        tol_NL = 1e-7
        maxiter_NL = 30

        FF_0 = VarForm_full(W,w,z,d,C_dam)     # Variational form tau_eps/u
        
        a_0 = lhs(FF_0)
        L_0 = rhs(FF_0)
        
        #Non linear loop
        while L2_error_u > tol_NL and it < maxiter_NL:
            it+=1
            #Compute solution 
            solve(a_0==L_0,w_,bcs,solver_parameters={'linear_solver': 'petsc',
                         'preconditioner': 'ilu'})
            
            #Damage propagation
            (u_inc,eps_inc) = w_.split(deepcopy=True) 
            dk.assign(d)  
            d_temp = calc_d(eps_inc,r,C0)
            d.assign(project(d_temp,Escalar))
            
            L2_error_u = assemble(((u_inc-u_old)**2)*dx)
            if it == maxiter_NL or np.isnan(L2_error_u): 
                exit()

            #Update variables
            r.assign(project(calc_new_r(eps_inc, r, C0), Escalar))    
            u_old.assign(u_inc)
            eps_old.assign(eps_inc) 
            
            #End non-linear loop
    
        #Project strains
        strain = project(eps_inc,ET2)
    
        # Save in pvd file
        d.rename('Damage index','d'); d_file << d
        u_inc.rename('Displacement','u'); u_file << u_inc
        strain.rename('Strain','Strain'); strain_file << strain  
        
        disp_x = u_inc((1,1))[0]  #max. displacement

        residual = action(a_0, w_)-L_0
        v_reac = Function(W)
        bc_Rx_i = DirichletBC(W.sub(0).sub(0),Constant(1.),contorno, 1)
        bc_Rx_i.apply(v_reac.vector())
        fx = assemble(action(residual, v_reac))  #Reaction X

        txt_file = open(file_path,'a')
        data = np.array([[disp_x,-fx]])
        np.savetxt(txt_file, data, delimiter=' ; ', newline='\n')
        txt_file.close()

        print('**************')
        print('Numero de step: ', inc)
        print('Error L2:', L2_error_u, 'Tolerancia' , tol_NL, 'Iteración', it, 'Máxima iteración' , maxiter_NL)
        print('**************')
        inc+=1

    return             
           
#%% MAIN
#import mesh
W,bcs,mesh,contorno = build_space()

#imposed displacement
load = Expression('t',degree=0,t=0.0)
bc_load = DirichletBC(W.sub(0).sub(0), load,contorno,2) #displacement X=1
bcs.append(bc_load)

# Scalar variables
E,nu = 28.3e3,0.18   # Young's Moduli [MPa], Poisson Coefficient
ft = 2.12            # tensile strength MPa
Gf = 0.373           # tensile fracture energy N/mm (69J/m^2)

L0 = 1               # Characteristic lenght
ce = 1               # Arbitrary constant for tau_eps
cu = 1               # Arbitrary constant for tau_u

dim = mesh.geometry().dim()  #model dimensions (2D -> 2)

L = E*Gf/(ft**2)         #Irwin's material length

# Creates folder to save results
d_file = File('test_mix/damage_index.pvd')
u_file = File('test_mix/displacement.pvd')
strain_file = File('test_mix/strain.pvd')

#Constant functions over the mesh
ET1 = TensorFunctionSpace(mesh, 'P', 1, shape=(3,3))
ET2 = TensorFunctionSpace(mesh, 'P', 1, shape=(2,2))
V0 = FunctionSpace(mesh,'P',1)
Escalar = FunctionSpace(mesh,'DP',0)     #Function space for scalar variables on the mesh
r0 = ft

#Verify if txt exist
file_path = 'data_mix.txt'
if os.path.isfile(file_path):
    os.remove(file_path)    #If txt exists -> remove

inc_f = 500
solve_damage(W,Escalar,bcs,mesh,contorno,r0,inc_f,file_path)

I know this is not actually an error, but the localization of the damage parameter doesn’t seems to be ok with the DP0projection. Do you think this would work better on Dolfinx?. I’m currently working on Dolfin v2019.1.0.

Regards!

1 Like
4

Sorry for the late reply. This code here is quite hard for me to parse, as I’m not familiar with damage models. Could you write out the actual strong formulation of your projection with LaTeX format?
Usually, that a DG-0 projection works and not a P-1 works, indicates that the function r is not in a P1 space. Are there any gradients in the formulation?

5

Well, the primal variables of my formulation (u,eps) are not the ones tha are bringing me problems, so I guess you are talking about the constitutive problem here

In that case, the damage model relates the stresses with the strains such that:

\sigma = (1-d)\mathbb{C}:\varepsilon

where d is an internal variable that its computed explicitly following:

d(r)=1-\frac{r_0}{r}\text{exp}\left( -2H_s \left( \frac{r-r_0}{r_0}\right) \right)

Here H_s and r_0 are material properties and r is another internal variable that evolves monotonically following:

r = \text{max}\{ r_0,\text{max}(\langle \overline{\sigma}_1\rangle) \}

and \langle \overline{\sigma}_1\rangle is the value of the largest principal effective stress (\overline\sigma=\mathbb{C}:\varepsilon) evaluated if the computed value is positive (\langle x \rangle=x, if x\ge0 and 0 otherwise).

I have included this in the nonlinear loop after the solver:

I guess the problem is that, even when d and r are not primal variables in the formulation, they are includded in the variational form (where \mathbb{C}(d)=(1-d)\mathbb{C}):

-(\gamma,\mathbb{C}(d):(\varepsilon-\nabla^s u))+\tau_\varepsilon(\gamma,\mathbb{C}(d):(\varepsilon-\nabla^s u)) - \tau_u(\text{div}(\mathbb{C}:\gamma),\text{div}(\mathbb{C}:\varepsilon)+f) = 0

(\nabla^s v,\mathbb{C}:\varepsilon)-\tau_\varepsilon(\nabla^s v,\mathbb{C}:(\varepsilon-\nabla^s u)) = (B,v) + (v,T)_{\partial\Omega}

Here the primal variables are u and \varepsilon. Also r depends on the primal variable \varepsilon that technically is a gradient (\varepsilon=\nabla^s u), but here it is computed as an independent variable at nodal level. Is there an alternative to include the values of d and r in the variational form without projecting?