Mesh Not updating in Adaptive refinement while solving coupled equations

9 months ago by

Hi, I am trying to solve coupled differential equations with staggered approach, for that I have tried to solved with adaptive refinement. Equation 1 is solved and mesh updated but Equation 2 is still solving o the previous mesh. I have used adapt(mesh.leaf_node()) to update the mesh. How can I resolve this issue.

Updated: Sorry I have not given full information (I thought full code will be lengthy).
Problem: equation 1 solved on updated mesh but equation 2 still solved on previous mesh why ?Figure 1: solution of equation 1 Figure 2: solution of equation 2

Thanks in advance

from dolfin import *
mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh,'CG',1)
W = VectorFunctionSpace(mesh,'CG',1)

p , q = TrialFunction(V), TestFunction(V)
u , v = TrialFunction(W), TestFunction(W)

l_fac = 2
Gc =  2.7
l_o = l_fac*mesh.hmin()
l_o = Expression('l_o', l_o=l_o)

lmbda = 121.15e3
mu = 85.77e3

def epsilon(v):
    return sym(grad(v))

def sigma(u):
    return 2.0*mu*epsilon(u) + lmbda*tr(epsilon(u))*Identity(2)

def en_dens(u):
    str_ele = 0.5*(grad(u) + grad(u).T)
    IC = tr(str_ele)
    ICC = tr(str_ele * str_ele)
    return (0.5*lmbda*IC**2) + mu*ICC

class top(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-10
        return abs(x[1]-1.0) < tol and on_boundary

class bottom(SubDomain):
    def inside(self,x,on_boundary):
        tol = 1e-10
        return abs(x[1]+0.0) < tol and on_boundary

Top = top()
Bottom = bottom()
u_Lx = Expression("t",t = 0.0)
bclx= DirichletBC(W.sub(0), Constant(0.0), Bottom)
bcly = DirichletBC(W.sub(1), Constant(0.0), Bottom)
bctx = DirichletBC(W.sub(1), u_Lx, Top)

bc_u = [bcly,bclx, bctx ]

u_old,u_conv, unew = Function(W),Function(W), Function(W)
alpha_old,alpha_conv = Function(V),Function(V)
# equation 1
E_du = ((1 - alpha_old)**2 + 1e-6)*inner(grad(v),sigma(u))*dx 
# equation 2
E_alpha = ( Gc*l_o*inner(grad(p),grad(q))+\
            ((Gc/l_o) + 2.*en_dens(unew))*inner(p,q)-\
u = Function(W)
p = Function(V)
ndim = mesh.geometry().dim()
L = inner(f, v)*dx 
problem_disp = LinearVariationalProblem(E_du,L,u,bc_u)

# error indicator
MF =  en_dens(u)*dx
solver_disp = AdaptiveLinearVariationalSolver(problem_disp, MF)
solver_disp.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
solver_disp.parameters["error_control"]["dual_variational_solver"]["symmetric"] = True

min_load = 0.0
max_load = 1.0

deltaT = 1e-3
t = deltaT
ut = 1

TOL = 1e-6

while t <= max_load:
    if t > 0.005:
        deltaT = 1e-4
    iter = 1
    toll = 1e-4
    err = 1
    maxiter = 100
    while err > toll:
        plot(mesh, key = 'mesh', title = 'updated mesh')
        hm = mesh.hmin()
        l_o = l_fac*hm
        print u.leaf_node().vector().array()
        err_u = errornorm(u.leaf_node(),u_old, norm_type = 'l2', mesh = None)
        err_alpha = errornorm(p,alpha_old, norm_type = 'l2', mesh = None)
        err = max(err_u,err_alpha)

        iter = iter + 1
        if err < toll:
            ux,uy = split(u.leaf_node())
            plot(ux,key = 'ux',title = 'u_dispx')
            plot(uy,key ='uy',title = 'u_dispy')
            plot(p.leaf_node(),range_min = 0.,range_max = 1.,key = 'alpha',title = 'alpha%.4f'%(t))

Community: FEniCS Project
You have not asked a question.  What is your question?
written 9 months ago by pf4d  
Mesh is not updated for the equation 2?
written 9 months ago by hirshikesh  
that is a statement with a question mark at the end...
written 9 months ago by pf4d  

1 Answer

9 months ago by
I don't know how you solve equation 2 since you didnt share that info with us.
However, maybe you need to adapt the problem for equation 2 with regard to the adapted mesh. Something like
adapt(Eu, mesh.leaf_node())
Eu = Eu.child()​
Thanks for the help, I will try this
written 9 months ago by hirshikesh  
I have tried but still mesh is not updating for the equation 2. I have given now updated MWE code (information: how i have solved equation 2).

Thanks for the help
written 8 months ago by hirshikesh  
You should be aware that you should use
plot(mesh.leaf_node(), key = 'mesh', title = 'updated mesh')​

to be sure to plot the adapted mesh. However, I think it is less confusing to use


so you don't have to deal with the leaf_node stuff at all.

Also there is propably still a

E_alpha = E_alpha.child()


written 8 months ago by Lukas O.  
Thanks for the help, after using this adapt(E_alpha,mesh.leaf_node())
I am getting Error: in method 'adapt', argument 1 of type 'dolfin::ErrorControl const &'
written 8 months ago by hirshikesh  
The following works for me
problem2 = LinearVariationalProblem(lhs(E_alpha),rhs(E_alpha),p)
solver2 = LinearVariationalSolver(problem2)
problem2 = adapt(problem2,mesh.leaf_node())
written 8 months ago by Lukas O.  
Thanks for the help. It works for me also ... :)
written 8 months ago by hirshikesh  
But Now concern is, Why mesh is refining everywhere (by half) ?
written 8 months ago by hirshikesh  
Because you did not specify any cell_markes, so it refines uniformly.
written 8 months ago by Lukas O.  
For the local refinement i have used error indicator MF in the code, I thought it will do the job. Is that something I am missing?
written 8 months ago by hirshikesh  
This is used to adapt the mesh internally in the adaptive solver. So you want to use the adapted mesh from the first solver also in the second one? I am not sure how to get the modified mesh from the adaptive solver. I think this is your real question.
written 8 months ago by Lukas O.  
yes, Thats what I wanted. Do we have any work around for this ?
written 8 months ago by hirshikesh  
or How can I define cell_markers based on the criteria en_dens() > some constant. I have tried this:
cell_markers = MeshFunction('bool', mesh)
for cell in cells(mesh):
    pp = cell.midpoint()
    if en_dens(pp) >= 10:

but getting error: Invalid type conversion: <class 'dolfin.cpp.mesh.Point'> can not be converted to any UFL type.

written 8 months ago by hirshikesh  
of course. en_dens takes a ufl function and you fed it with a point.
written 8 months ago by Lukas O.  
If you just delete adapt(mesh.leaf_node()) you should be fine.
written 8 months ago by Lukas O.  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »