Integrating on subdomain

Dear fellows,
I am trying to solve the finite element problem on the marked domain while keeping the solution intact on the unmarked domain. I have marked the concerned domain with “mark_diff_domain” as domain_id of “1”. The code runs fine if I do not specify marked subdomain and integrate over complete domain using simply “dxx”.
My objective is to integrate over marked subdomain only and keep the unmarked domain’s solution intact until future time steps, when “mark_diff_domain” will be updated to include some other cells from currently unmarked region into solution domain.
So, now when I try to integrate and solve my variational form only on this domain using dxx(1) then it retrieves the following error:

sudo python3 thermal_analysis_mwe.py mesh_to_import.xml
step: 1 0.01
Traceback (most recent call last):
File “thermal_analysis_mwe.py”, line 86, in
solve(F == 0, u, bcs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 220, in solve
_solve_varproblem(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py”, line 266, in _solve_varproblem
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at


*** fenics-support@googlegroups.com


*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.


*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0


*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Can someone please suggest me the solution or correct me in my syntax if I am missing something important? I attach the syntax file in attachment. Thank you in advance!

Link to the code and mesh file:

https://drive.google.com/file/d/1fB4146_klPikr1RsDE1XWzsIP-pDg8tu/view?usp=sharing,%20https://drive.google.com/file/d/1gO_DxWujhZ02NkSSJdAyBtAT8jMay1Sd/view?usp=sharing

https://drive.google.com/file/d/1gO_DxWujhZ02NkSSJdAyBtAT8jMay1Sd/view?usp=sharing

Dear mashsquad777

It would be better if you put the code herein, since your link is not working. Generally speaking, the whole domain needs to be described, otherwise, it will not converge. But there is also

cells_mat = MeshFunction(‘size_t’, mesh, 3, 0)
mesh_markedbyone = MeshView.create(cells_mat, 1)

so you can define on this new mesh an integral measure and solve a problem.

Best, Emek

Dear Emek,
Thank you very much for your guiding response. Is "MeshView.create(cells_mat,1) is same as the Submesh feature? So the idea is to creates mesh of subdomain and that subdomain can be used for solution by integrating on whole subdomain mesh instead of integrating on partial subdomain e.g. dxx(1)?.
Sorry if the link is not accessible. Following is the code I am working upon:

from fenics import *
import logging
set_log_level(logging.WARNING)
import sys


#'Unit cube mesh

mesh_final = UnitCubeMesh(8,8,8);



# Define boundary subdomains
tol_boundary = 1e-14

class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0, tol_boundary)
class BoundaryY1_diff(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1, tol_boundary)
class DiffusionHalfDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.5))
        

# Time-stepping
t = 0
istep = 0
dt = 1e-2
num_steps = 250

# Functional space
V = FunctionSpace(mesh_final, 'P', 1)

# Define initial value

u_n = interpolate(Constant(400), V)
u = Function(V)
u.interpolate(Constant(400))
v = TestFunction(V)
f = Constant(0)

# marking the subdomain 
mark_diff_domain = MeshFunction("size_t", mesh_final, mesh_final.topology().dim(), 0)
mark_domain = DiffusionHalfDomain()
mark_diff_domain.set_all(0)
mark_domain.mark(mark_diff_domain, 1)


for n in range(num_steps):

    # Update current time
    t += dt

    istep += 1
    print("step:", istep, t)

    # Mark boundaries 
    boundary_markers = MeshFunction("size_t", mesh_final, mesh_final.topology().dim()-1, 9999)
    by0 = BoundaryY0(); by0.mark(boundary_markers, 2)
    by1 = BoundaryY1_diff(); by1.mark(boundary_markers, 3)

    #Dirichlet conditions
    bcs = []

    bc_cur = DirichletBC(V, Constant(400), boundary_markers, 2)
    bcs.append(bc_cur)
    bc_cur = DirichletBC(V, Constant(500), boundary_markers, 3)
    bcs.append(bc_cur)
    
    # physical parameters
    cpv = 50
    rho = 4420
    cpvrho = cpv*rho
    tdiff = 1000

    #dx = Measure('dx', domain = mesh_final, subdomain_data = mark_diff_domain, subdomain_id = 1)
    dx = Measure("dx")
    dxx = dx(subdomain_data=mark_diff_domain)
    
    F = cpvrho*u*v*dxx(1) + tdiff*dt*dot(grad(u), grad(v))*dxx(1) - cpvrho*(u_n + dt*f)*v*dxx(1)       

    # Compute solution
    solve(F == 0, u, bcs)

    # previous solution
    u_n.assign(u)

    # Plot solution
    xdmff = XDMFFile("forttemp" + format(istep, '04') + ".xdmf")
    u.rename("u","temp")
    xdmff.write(u)
    xdmff.close()

Please format your code using ``` encapsulation, and make sure all indentation is correct, such that anyone can copy-paste your code and run it without modifications.

Sorry @dokken I am new and did not know about this feature. However I encapsulated the code now. Thanks!

Yes, it is practically the same as submesh, but it is also working in parallel computing very nicely. It is mostly used for projecting the solution on the whole (parent) mesh to a submesh (MeshView.create) and show it nicely in ParaView. But you can also use it for solving a problem, if the measure is built by that mesh only! You cannot solve on a part of the mesh like dxx(1) it should be a new measure and function space on this new submesh. But after solving it, you can use project to exchange information.

Best, Emek

Thanks a lot Emek! It solves now part of my problem here. :slight_smile:
Other than thermal analysis, I am also using fenics solid mechanics application and my subdomain is evolving at every time or so called load step.
So if I use submesh there to make new updated subdomain and functional space from it every time then the problem gets redefined and unfortunately the fsm nonlinear solver then forgets about the history variables in old domain and solves for the current (updated submesh).
It would be great if I can keep the history variables like plastic strain, displacements etc. from the previous subdomain to the updated subdomain zones where it was solved in previous step.
I hope I described the objective clearly.

Not necessarily, it depends on the implementation and if you update the solution by using assign(…) it is not renewing the whole assembly.

If you can have a minimal working example, where the problem happens, post it here.

Best, Emek

Thanks Emek for the idea of posting mwe and also your on going support.
So far in current simple code the analysis resets on each loop iteration when the domain evolves and subdomain is updated.
So my objective is that the fsm should remember what happened in the current solution subdomain e.g. history variables of the plastic strains and deformations etc. so that when next solution domain is activated e.g. quarter subdomain then the material behaves as it is already loaded when only half domain was active. Here it is very simple example but in my actual application even the loading on previously activated subdomain will also vary in future loading steps and certain parts of domain will have to behave in a certain way depending if they were already plasticaly deformed or not.

Following is the simple testing code from fsm application where I am trying to sequentially evolve the domain of unit cube in x direction from half, quarter to whole cubic domain and solving with prescribed displacement boundary condition on top (y=1):

from dolfin import *
import fsm

#set_log_level(10)

class DirichletBoundaryX(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] < DOLFIN_EPS)
class DirichletBoundaryY(SubDomain):
    def inside(self, x, on_boundary):
        return (x[1] < DOLFIN_EPS)
class DirichletBoundaryZ(SubDomain):
    def inside(self, x, on_boundary):
        return (x[2] < DOLFIN_EPS)

class PrescribedDisplacementY(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (1-DOLFIN_EPS, 1) ) 

class HalfDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.5))
class QuarterDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.75))
class WholeDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 1))

# making output file
file1 = XDMFFile('output.xdmf')      

# Making complete mesh
mesh_Complete = UnitCubeMesh(8,8,8);
mark_sub_domain = MeshFunction("size_t", mesh_Complete, mesh_Complete.topology().dim(), 0)

number_of_submesh = 3
for submesh_number in range(number_of_submesh): 
    
    if (submesh_number == 0) :
        mark_domain = HalfDomain()
        print('HalfDomain')
    elif (submesh_number == 1) :
        mark_domain = QuarterDomain()
        print('QuarterDomain')
    elif (submesh_number == 2) :
        mark_domain = WholeDomain()
        print('WholeDomain')
    
    #mark_domain = WholeDomain()
    # marking the subdomain 
    
    mark_sub_domain.set_all(0)
    mark_domain.mark(mark_sub_domain, 1)

    mesh = SubMesh(mesh_Complete, mark_sub_domain, 1)
    
    E = 20000.0;
    nu = 0.3;

    scheme = "default"
    degree = 3
    dx = Measure("dx")
    dx = dx(degree=degree, scheme=scheme)

    V  = VectorFunctionSpace(mesh, "Lagrange", 2)
    element_t = VectorElement("Quadrature", mesh.ufl_cell(), degree=3, dim=36, quad_scheme=scheme)
    Vt = FunctionSpace(mesh, element_t)
    element_s = VectorElement("Quadrature", mesh.ufl_cell(), degree=3, dim=6, quad_scheme=scheme)
    Vs = FunctionSpace(mesh, element_s)
    
    zero = Constant(0.0)
    prescribed_displacement = Constant(1e-3*(submesh_number+1))

    bc0 = DirichletBC(V.sub(0), zero, DirichletBoundaryY(), method="pointwise")
    bc1 = DirichletBC(V.sub(1), zero, DirichletBoundaryY(), method="pointwise")
    bc2 = DirichletBC(V.sub(2), zero, DirichletBoundaryY(), method="pointwise")
    bc3 = DirichletBC(V.sub(1), prescribed_displacement, PrescribedDisplacementY(), method="pointwise")

    bcs = [bc0, bc1, bc2, bc3]

    E_t = 0.3*E
    hardening_parameter = E_t/(1.0 - E_t/E)
    yield_stress = 9.0

    u = Function(V, name="u")

    def eps(u):
        return as_vector([u[i].dx(i) for i in range(3)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1), (0, 2), (1, 2)]])

    def sigma(s):
        #s = ss.function_space()
        return as_matrix([[s[0], s[3], s[4]], [s[3], s[1], s[5]], [s[4], s[5], s[2]]])

    def tangent(t):
        #t = tt.function_space()
        return as_matrix([[t[i*6 + j] for j in range(6)] for i in range(6)])

    J2 = fsm.python.cpp.plasticity_model.VonMises(E, nu, yield_stress, hardening_parameter)
    Qdef = fsm.UFLQuadratureFunction(eps(u), element_s, mesh)
    fsm_constitutive_update = fsm.ConstitutiveUpdate(Qdef, J2)
    #fsm_tangent = QuadratureFunction(mesh, Vt.element(), fsm_constitutive_update, fsm_constitutive_update.w_tangent())
    #fsm_stress = QuadratureFunction(mesh, Vs.element(), fsm_constitutive_update.w_stress())
    fsm_tangent = fsm.QuadratureFunction(Vt, fsm_constitutive_update.w_tangent(), fsm_constitutive_update)
    fsm_stress  = fsm.QuadratureFunction(Vs, fsm_constitutive_update.w_stress())


    v = TestFunction(V)
    uTrial = TrialFunction(V)

    a = inner(eps(v), dot(tangent(fsm_tangent), eps(uTrial)) )*dx
    L = inner(grad(v), sigma(fsm_stress))*dx

    nonlinear_problem = fsm.PlasticityProblem(a, L, u, fsm_tangent, fsm_stress, bcs)

    nonlinear_solver = NewtonSolver()
    nonlinear_solver.parameters["convergence_criterion"] = "incremental";
    nonlinear_solver.parameters["maximum_iterations"]    = 50;
    nonlinear_solver.parameters["relative_tolerance"]    = 1.0e-6;
    nonlinear_solver.parameters["absolute_tolerance"]    = 1.0e-15;

    # File names for output

    eps_p_eq = fsm_constitutive_update.eps_p_eq()
    #fsm_constitutive_update.eps_p_eq().compute_mean(eps_eq);

    element_eps_p_eq_project = FiniteElement("CG", mesh.ufl_cell(), degree=1)
    V_eps_p_eq_project = FunctionSpace(mesh, element_eps_p_eq_project)
    eps_p_eq_project = Function(V_eps_p_eq_project, name="eps_p_eq")

    # Solve non-linear problem
    nonlinear_solver.solve(nonlinear_problem.cpp_object(), u.vector());

    # Update variables
    fsm_constitutive_update.update_history();

    # Write output to files
    file1.write(u, t=float(submesh_number));
    eps_p_eq_project.vector()[:] = project(eps_p_eq, V_eps_p_eq_project, solver_type='gmres',
                                        form_compiler_parameters={
                                            "representation": parameters["form_compiler"]["representation"],
                                            "quadrature_scheme": scheme,
                                            "quadrature_degree": degree
                                        }
                                    ).vector()
    file1.write(eps_p_eq_project, t=float(submesh_number));
file1.close()

Well, technically it is not a minimal working example since you use non-standard libraries like fsm. I may guess that your

eps_p_eq_project = Function(V_eps_p_eq_project, name="eps_p_eq")
eps_p_eq_project.vector()[:] = project(...)

is causing a recompilation. Try to define

eps_p_eq_project = Function(V_eps_p_eq_project, name="eps_p_eq")

outside the loop and then within the loop use

assign(eps_p_eq_project, project(...) )

In this way you make sure that you use the same object and only exchange its numerical values. So the FFC does not recompile, since all objects are still the same.

Best, Emek

Thanks a lot Emek. I tried as you suggested. It is working fine for keeping the history result of previous time or loading steps. Now only thing remaining is that updating the submesh as the domain evolves. As far as my limited experience goes, if I assign the updated mesh to the fenics problem definition then it may again reset the system for each updated mesh.


from dolfin import *
import fsm

#set_log_level(10)

class DirichletBoundaryX(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] < DOLFIN_EPS)
class DirichletBoundaryY(SubDomain):
    def inside(self, x, on_boundary):
        return (x[1] < DOLFIN_EPS)
class DirichletBoundaryZ(SubDomain):
    def inside(self, x, on_boundary):
        return (x[2] < DOLFIN_EPS)

class PrescribedDisplacementY(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[1], (1-DOLFIN_EPS, 1) ) 

class HalfDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.5))
class QuarterDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 0.75))
class WholeDomain(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0, 1))

# making output file
file1 = XDMFFile('output.xdmf')      
file1.parameters['functions_share_mesh'] = True
file1.parameters['rewrite_function_mesh'] = False
file1.parameters["flush_output"] = True

# Making complete mesh
mesh_Complete = UnitCubeMesh(8,8,8);
mark_sub_domain = MeshFunction("size_t", mesh_Complete, mesh_Complete.topology().dim(), 0)

number_of_submesh = 3
number_of_steps = 3
# marking the subdomain
    
mesh = mesh_Complete 

E = 20000.0;
nu = 0.3;

scheme = "default"
degree = 3
dx = Measure("dx")
dx = dx(degree=degree, scheme=scheme)

V  = VectorFunctionSpace(mesh, "Lagrange", 2)
element_t = VectorElement("Quadrature", mesh.ufl_cell(), degree=3, dim=36, quad_scheme=scheme)
Vt = FunctionSpace(mesh, element_t)
element_s = VectorElement("Quadrature", mesh.ufl_cell(), degree=3, dim=6, quad_scheme=scheme)
Vs = FunctionSpace(mesh, element_s)

zero = Constant(0.0)
prescribed_displacement = Expression(" 5e-3*t1",t1=0,degree=2)

bc0 = DirichletBC(V.sub(0), zero, DirichletBoundaryY(), method="pointwise")
bc1 = DirichletBC(V.sub(1), zero, DirichletBoundaryY(), method="pointwise")
bc2 = DirichletBC(V.sub(2), zero, DirichletBoundaryY(), method="pointwise")
bc3 = DirichletBC(V.sub(1), prescribed_displacement, PrescribedDisplacementY(), method="pointwise")

bcs = [bc0, bc1, bc2, bc3]

E_t = 0.3*E
hardening_parameter = E_t/(1.0 - E_t/E)
yield_stress = 9.0

u = Function(V, name="u")

def eps(u):
    return as_vector([u[i].dx(i) for i in range(3)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1), (0, 2), (1, 2)]])

def sigma(s):
    #s = ss.function_space()
    return as_matrix([[s[0], s[3], s[4]], [s[3], s[1], s[5]], [s[4], s[5], s[2]]])

def tangent(t):
    #t = tt.function_space()
    return as_matrix([[t[i*6 + j] for j in range(6)] for i in range(6)])

J2 = fsm.python.cpp.plasticity_model.VonMises(E, nu, yield_stress, hardening_parameter)
Qdef = fsm.UFLQuadratureFunction(eps(u), element_s, mesh)
fsm_constitutive_update = fsm.ConstitutiveUpdate(Qdef, J2)
#fsm_tangent = QuadratureFunction(mesh, Vt.element(), fsm_constitutive_update, fsm_constitutive_update.w_tangent())
#fsm_stress = QuadratureFunction(mesh, Vs.element(), fsm_constitutive_update.w_stress())
fsm_tangent = fsm.QuadratureFunction(Vt, fsm_constitutive_update.w_tangent(), fsm_constitutive_update)
fsm_stress  = fsm.QuadratureFunction(Vs, fsm_constitutive_update.w_stress())


v = TestFunction(V)
uTrial = TrialFunction(V)

a = inner(eps(v), dot(tangent(fsm_tangent), eps(uTrial)) )*dx
L = inner(grad(v), sigma(fsm_stress))*dx

nonlinear_problem = fsm.PlasticityProblem(a, L, u, fsm_tangent, fsm_stress, bcs)

nonlinear_solver = NewtonSolver()
nonlinear_solver.parameters["convergence_criterion"] = "incremental";
nonlinear_solver.parameters["maximum_iterations"]    = 50;
nonlinear_solver.parameters["relative_tolerance"]    = 1.0e-6;
nonlinear_solver.parameters["absolute_tolerance"]    = 1.0e-15;

# File names for output

eps_p_eq = fsm_constitutive_update.eps_p_eq()
#fsm_constitutive_update.eps_p_eq().compute_mean(eps_eq);

element_eps_p_eq_project = FiniteElement("CG", mesh.ufl_cell(), degree=1)
V_eps_p_eq_project = FunctionSpace(mesh, element_eps_p_eq_project)
eps_p_eq_project = Function(V_eps_p_eq_project, name="eps_p_eq")

for step_number in range(number_of_steps): 
    
    #assign(mesh, mesh_collection[step_number])
    
    if (step_number == 1):
        prescribed_displacement.t1 = 1
    else:
        prescribed_displacement.t1 = 0

    # Solve non-linear problem
    nonlinear_solver.solve(nonlinear_problem.cpp_object(), u.vector());

    # Update variables
    fsm_constitutive_update.update_history();

    # Write output to files
    file1.write(u, t=float(step_number));

    assign( eps_p_eq_project, project(eps_p_eq, V_eps_p_eq_project, solver_type='gmres',
                                        form_compiler_parameters={
                                            "representation": parameters["form_compiler"]["representation"],
                                            "quadrature_scheme": scheme,
                                            "quadrature_degree": degree
                                        }
                                    ) )
                                    
    file1.write(eps_p_eq_project, t=float(step_number));
    

file1.close()

Yes, if you change the topology, the matrices need to be re-calculated. But if you change the coordinates of the mesh, as we do in Fluid-Structure interaction (moving mesh), then the topology remains the same and no re-assembly takes place.

Thanks for sharing the information. I will keep that in mind. Then in my case as the domain evolves, so I think I have to project the old solution on new subdomain at common zones declaring uncommon zones (newly introduced in result of evolving of domain) solution as zero.
Is there a way in your knowledge to do so or maybe another trick? As you know I am using fsm based code.
So far I also tried using transfer_matrix = PETScDMCollection.create_transfer_matrix(V_old_domain,V_new_domain) for displacement solution projection in new submesh before handing problem over to nonlinear solver but this projects whole solution to new subdomain even to newly evolved zones where solution is extrapolated as non zero.