Preserve Material discontinuities in projection for postprocessing

1

Hi,

I am using Fenics to perform Finite-Element simulations for coupled magneto-mechanical problems (stationary) and have a problem with the calculation of postprocessing quantities for visualization in Paraview. In my case, the domain contains subdomains with different material parameters. Assigning those parameters and solving the problem is no problem, but when it comes to postprocessing, I always get smoothed results at the material discontinuities which I would like to prevent.

I thought using DG spaces for the projection of my postprocessing quantities would help… Indeed, it does: as long as I am using a DG0 space. However, since my solution for the primary displacement and magnetic potential fields is second order, I would like to be able to also use linear function spaces for all quantities that involve gradients, etc. If I try to do a projection (local or global - does not change the effect) with DG1 elements, I end up with results that are smoothed over the discontinuity which is not surprising but definitely not what I want!

Is there any possibility to generate postprocessing results for Paraview files in which the discontinuities between different subdomains are preserved?

Here is what I basically do:

class PostProcessingQuantity(df.Function):

    #########################
    # Initialization method #
    #########################
    def __init__(self,function,space,name,file):

        # initialize parent Function object
        super().__init__(space,name=name)

        # save required variables
        self.W = space                                                          # Function space for projections
        self.dx = ufl.dx(self.W.mesh())                                         # measure dx (required for local projections)
        self.u = function                                                       # function to project
        self.file=file

        # distinguish continuous and discontinuous projection spaces
        if self.W.ufl_element().family() == "Discontinuous Lagrange":           # projection space is discontinuous
            self.localProjection=True                                           # -> enable local projection
            self.du = df.TrialFunction(self.W)                                  # -> initialize TrialFunction du
            self.v = df.TestFunction(self.W)                                    # -> initialize TestFunction v
        else:                                                                   # projection space is continuous
            self.localProjection=False                                          # -> do not enable local projection


    def update(self,time):
        self.performProjection()
        self.writeToResultFile(time)



    #######################################################
    # Method for projection to the defined solution space #
    #######################################################
    def performProjection(self):

        # distinguish global and local projections
        if self.localProjection:                                                # local projection (discontinuous projection spaces)
            lhs = df.inner(self.du,self.v)*self.dx                              # -> determine left-hand-side of system to solve
            rhs = df.inner(self.u,self.v)*self.dx                               # -> determine right-hand-side of system to solve
            solver = df.LocalSolver(lhs,rhs)                                    # -> define solver to be local, i.e. element-wise
            solver.factorize()                                                  # -> use factorization to prevent unnecessary calculations
            solver.solve_local_rhs(self)                                        # -> solve local system and assign solution to PostProcessingQuantity
        else:                                                                   # global projection
            self.assign(df.project(self.u,self.W))                              # -> solve global system and assign solution to PostProcessingQuantity


    def writeToResultFile(self,time):
        self.file.write(self,time)

Performing the projection on a per-subdomain basis also does not change the effect of averaged quantities over the material interface. Maybe, I shouldn’t be using projection at all, but I was not able to figure out, how I could write a quantity like grad(u) - where u is one of my primary fields for which I solve the problem - to a Paraview file without projecting

2

Are you using XDMF and write_checkpoint instead of write function ?

1 Like
3

To add some context, see for instance: Access function in time dependent XDMFFile

1 Like
4

I tried XDMF and pvd files. The results did not change… But I never tried to use the write_checkpoint function, only the write function. I will give it a try and see, if that already answers my question

5

OK, this did the trick. Now the solutions look more reasonable.

As a follow-up question: I still have problems writing my primary fields to an xdmf-file. Currently I am using periodic boundary conditions to determine the fluctuations for my primary field which is approximated by CG2 ansatz functions. When I want to save the complete displacement field for post-processing, I have to add the average displacement which is prescribed via an expression. So what I do is basically the following:

dot((FMacro-I),x)+u

Here, FMacro is an Expression and x is is defined as SpatialCoordinate. With this Expression, I cannot simply use interpolation to write the total displacement to an XDMFFile - this does not work since the equation above has no cpp_object attribute. However, if I use projection to project the solution to the same space my primary fluctuation field u lives in, the results are complete nonsense. I get diverging results as if the projection problem was ill-posed.

Is there any other hint on how I can accomplish the post-processing of my primary fields within a solution space that is the same as my finite element simulations (DG1 does work, CG1 gives a spotted solution that apparently also is not what I want, see the attached image)?

result
result1922Ă—1175 760 KB

Except using the write_checkpoint method of XDMFFiles, the code for postprocessing has not changed…

6

So what happens if you only try to project dot(FMacro-I, x) and visualize it with write_checkpoint? As you have not given a minimal example illustrating what FMacro is, there isn’t much I can do.
Also, what happens if you only project u into its own space? Does that work?

7

Sorry, I should have uploaded more of the code… Yes, I checked projecting only u or dot(FMacro-I,x) and both give similar results to what is shown in the picture above. If I just use u without projecting and write it to a .pvd-file, the solution looks as expected. It’s a little bit tricky to really give a minimal working example, since I have written most of the codes features to separate python classes to be able to reuse them for different problems, but I will give it a try:

##############################################################################
# TEST FILE FOR THE IMPLEMENTATION OF PERIODIC BOUNDARY CONDITIONS IN FENICS #
##############################################################################

# load required libraries
#########################
from dolfin import *
import numpy as np

# Optimization options for the form compiler
############################################
parameters["form_compiler"]["cpp_optimize"] = True
parameters ["form_compiler"]["quadrature_degree"] = 2
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


#######################################
# Classes with reusable code features #
#######################################
class PeriodicDomain(SubDomain):

    def __init__(self,mesh,boundingBox,periodicAxes=None,tolerance=DOLFIN_EPS):

        # plausibility checks
        if mesh is None:
            raise TypeError("No mesh object for the definition of periodicity constraints passed. Check your input data.")
        if boundingBox is None:
            raise TypeError("No RVE bounding box (minimum and maximum points) for the definition of periodicity constraints passed. Check your input data.")
        elif any(np.array(boundingBox[0]) > np.array(boundingBox[1])):
            raise ValueError("Impossible mesh bounding box passed for the definition of periodicity constraints. In the mesh bounding box, no coordinate of the minimum point can be larger than the corresponding maximum points coordinate. Check your input data.")

        # calculate relevant information from bounding box array
        axesMin=np.array(boundingBox[0])
        axesMax=np.array(boundingBox[1])
        relevantAxes=np.flatnonzero(axesMax-axesMin)
        dimension=len(relevantAxes)
        center=Point((axesMax[relevantAxes]+axesMin[relevantAxes])/2)

        # plausibility checks for periodicAxes
        if periodicAxes is None:                                                # check if periodicAxes has a value
            raise TypeError("Variable \"periodicAxes\" not set for calculation of periodicity constraints! Check your input data.")
        elif len(np.shape(periodicAxes)) > 1:                                   # check for right amount of dimensions
            raise ValueError("Wrong amount of array dimensions for variable \"periodicAxes\"! For a cuboid periodic domain, the variable \"periodicAxes\" can only be one-dimensional. Check your input data.")
        elif np.any(~np.isin(periodicAxes,relevantAxes)):
            raise ValueError("Detected intended periodicity constraint for an axis that is either outside [0,1,2], i.e. not a spatial axis, or an axis in which the mesh has zero extent. Only spatial axes in which the mesh has a nonzero extent can be set periodic. Check your input data.")

        # initialize parent class
        super().__init__()

        # set relevant information for further calculations
        self.fixedPoint=PointDomain(mesh,center,tolerance)                      # get subdomain for fixed center point (closest vertex)
        self.relevantAxes=relevantAxes                                          # relevant calculation axes
        self.periodicAxes=periodicAxes                                          # periodic axes
        self.axesMin=axesMin                                                    # minimum axes coordinates (lower, left, back point of RVE)
        self.axesMax=axesMax                                                    # maximum axes coordinates (upper, right, front point of RVE)
        self.tolerance=tolerance                                                # allowed tolerance


    ##########################################################################
    # Overwrite inside method of SubDomain class to define master subdomains #
    ##########################################################################
    def inside(self,x,on_boundary):
        """
        This method defines the master subdomains which here are the negative
        boundary entities without the parts that also belong to a positive
        boundary entity. This allows to apply constraints without any duplicate
        information. Nodes that are not part of a master entity but also do not
        belong to a slave entity - since, e.g., only selected directions are
        supposed to be periodic - are handled within the mapping method.
        """
        # check if point x belongs to AT LEAST ONE master but NO slave entities
        xInMasterEntity=False                                                   # initialze xInMasterEntity flag to be False by default
        for axis in self.periodicAxes:                                          # loop over all periodic axes
            if abs(x[axis]-self.axesMax[axis]) < self.tolerance:                # -> x is part of the current axis slave entity
                xInMasterEntity=False
                break                                                           # ->-> stop further checks: x cannot be part of a master entity
            elif abs(x[axis]-self.axesMin[axis]) < self.tolerance:              # -> x is part of the current axis master entity
                xInMasterEntity=True                                            # ->-> set master entity membership information

        return xInMasterEntity


    #################################################################################
    # Overwrite map method of SubDomain class to define the mapping master -> slave #
    #################################################################################
    def map(self,x,y):
        """
        This method defines how the mapping between nodes on the slave and master
        entities is performed. Here, coordinates x that are identified to be part
        of a slave domain are mapped to the master domain with coordinates y.
        Points that do neither belong to a master nor to a slave entity, are
        handled with the simple mapping y=x, i.e. their coordinates remain
        unchanged.
        """
        xInSlaveEntity=False                                                    # initialize xInSlaveEntity flag to be False by default
        for axis in self.relevantAxes:                                          # loop over all axes with nonzero extent
            if axis in self.periodicAxes:                                       # current axis is supposed to be periodic
                if abs(x[axis]-self.axesMax[axis]) < self.tolerance:            # -> point is on slave entity for the current axis
                    y[axis]=self.axesMin[axis]                                  # ->-> set coordinate to master nodes coordinate
                    xInSlaveEntity=True                                         # ->-> update slave entity membership flag
                else:                                                           # -> point is not on slave entity for the current axis
                    y[axis]=x[axis]                                             # ->-> leave coordinate unchanged
            else:                                                               # current axis is not supposed to be periodic
                y[axis]=x[axis]                                                 # -> leave coordinate unchanged



#####################
# PointDomain class #
#####################
# This class implements a single point domain as a child of the FEniCS-internal
# subdomain class: the parent classes inside method is overwritten to find the
# closest vertex to the specified point.
class PointDomain(SubDomain):

    def __init__(self,mesh,point,tolerance=DOLFIN_EPS):

        # plausibility checks
        if mesh is None:
            raise TypeError("No mesh object passed for PointDomain passed. Check your input.")
        if point is None:
            raise TypeError("No point for the definition of single point subdomain passed. Check your input data.")

        # get mesh bounding box tree
        bbt=mesh.bounding_box_tree()

        # initialize parent SubDomain class
        super().__init__()

        # check if point is within the current mesh
        if len(bbt.compute_collisions(point)) > 0:                              # point is in this part of the mesh
            closestCell,_=bbt.compute_closest_entity(point)                     # get closest entity to point, i.e. cell that includes the point
            minimumDistance=1/tolerance                                         # initialize arbitrary big number
            for vertex in vertices(Cell(mesh,closestCell)):                     # loop over vertices of the closest cell
                distanceToPoint=point.distance(vertex.point())                  # -> get distance of vertex to cell
                if distanceToPoint < minimumDistance:                           # -> check if distance is below minimum distance
                    closestVertex=vertex.point()                                # ->-> save new closest vertex
                    minimumDistance=distanceToPoint                             # ->-> save new minimum distance
        else:                                                                   # point is not in this part of the mesh
            closestVertex=None                                                  # assign None to closest vertex

        # save relevant information
        self.tolerance=tolerance
        self.closestVertex=closestVertex


    ###################################################################################
    # Overwrite inside method of SubDomain class to identify nodes inside PointDomain #
    ###################################################################################
    def inside(self,x,on_boundary):
        if self.closestVertex is not None:                                      # closest vertex is not None
            # only return True if all coordinates are close to the specified vertex
            for axis in range(0,len(x)):
                if not abs(x[axis]-self.closestVertex.array()[axis]) < self.tolerance:
                    return False
            return True
        else:                                                                   # closest vertex is None
            return False



################################################################################
# >>>>>>> PROBLEM  TO SOLVE STARTS HERE
################################################################################

# relevant simulation parameters:
#################################
# domain information
domainBB=[[0, 0, 0],[1, 1, 0]]                                                  # domain bounding box

# time
nIncs=5
dt=1/nIncs

# loading
d11=Constant(0.5)

# elastic parameters
E = 2e6
nu = 0.48

# Lame constants
mu = Constant(E/(2*(1+nu)))
lam = Constant(E*nu/((1+nu)*(1-2*nu)))

# domain information (mesh, function spaces, ...)
#################################################
# import mesh (replaced with built-in meshing functionality for testing purposes)
mesh = UnitSquareMesh(10, 10)
dx=Measure("dx",domain=mesh)
x=SpatialCoordinate(mesh)

# generate subdomains for periodicity constraints and fixed center point
periodicDomain=PeriodicDomain(mesh,domainBB,periodicAxes=[0,1],tolerance=1e-3)  # periodicity constraints

# define function spaces
V = VectorElement("P", mesh.ufl_cell(), 2)
W = FunctionSpace(mesh,V,constrained_domain=periodicDomain)

# define relevant functions
u=Function(W)
w_u=TestFunction(W)
du=TrialFunction(W)

# define boundary conditions (fixed point)
BC=DirichletBC(W,Constant((0.,0.)),periodicDomain.fixedPoint,method="pointwise")

# required mathematical expressions
###################################
FMacro = Expression((("1+0.5*t", "0"),("0", "1/(1+0.5*t)")), t=0, degree=3 )

# Kinematics
############
dim = u.geometric_dimension()
I = Identity(dim)
F = FMacro + grad(u)
F1 = inv(F)
b = F*F.T
detF = det(F)

# Kinetics
##########
P = F1*( (mu*(b-I) + lam/2*(detF*detF-1)*I)) # standard Neo-Hookean approach

# Define weak form
##################
WF = inner(P.T,grad(w_u))*dx
J = derivative(WF,u,du)

# Setup solver
##############
problem = NonlinearVariationalProblem(WF, u, BC, J=J,
                                      form_compiler_parameters=ffc_options)
solver  = NonlinearVariationalSolver(problem)
solverParams=solver.parameters
solverParams['newton_solver']['absolute_tolerance'] = 1E-5
solverParams['newton_solver']['relative_tolerance'] = 1E-4
solverParams['newton_solver']['maximum_iterations'] = 20
solverParams['newton_solver']['relaxation_parameter'] = 1.0
solverParams['newton_solver']['linear_solver'] = 'mumps'


# prepare post-processing
#########################
# spaces
WPost = VectorFunctionSpace(mesh,"CG",2)

# functions
uAvg = Function(WPost,name="AverageDisplacement")
uFluc = Function(WPost, name="FluctuationDisplacement")
uTot = Function(WPost, name="Displacement")

# files without context manager
fileFluc=File("Solution/Fluctuation.pvd")

# Incremental solution
######################
t=0

# initial result generation for XDMF Files
uAvg.assign(project(dot(FMacro-I,x),V=WPost))
uFluc.assign(project(u,V=WPost))
uTot.assign(project(dot(FMacro-I,x)+u,V=WPost))
with XDMFFile("Solution/AverageDisplacement.xdmf") as file:
    file.write_checkpoint(uAvg,"AverageDisplacement",t)
with XDMFFile("Solution/FluctuationDisplacement.xdmf") as file:
    file.write_checkpoint(uFluc,"FluctuationDisplacement",t)
with XDMFFile("Solution/Displacement.xdmf") as file:
    file.write_checkpoint(uTot,"Displacement",t)

for i in range(nIncs):

    # Update time and BCs
    t += dt
    FMacro.t = t

    # solution
    solver.solve()

    # projection
    uAvg.assign(project(dot(FMacro-I,x),V=WPost))
    uFluc.assign(project(u,V=WPost))
    uTot.assign(project(dot(FMacro-I,x)+u,V=WPost))

    # XDMFFile generation
    with XDMFFile("Solution/AverageDisplacement.xdmf") as file:
        file.write_checkpoint(uAvg,"AverageDisplacement",t,file.Encoding.HDF5,append=True)
    with XDMFFile("Solution/FluctuationDisplacement.xdmf") as file:
        file.write_checkpoint(uFluc,"FluctuationDisplacement",t,file.Encoding.HDF5,append=True)
    with XDMFFile("Solution/Displacement.xdmf") as file:
        file.write_checkpoint(uTot,"Displacement",t,file.Encoding.HDF5,append=True)

This should capture the behavior I am observing. If I use a DG1 space for WPost, everything works just fine, with the applied CG2 space, the results are as described above…

8

The quadrature degree is too low. Removing this parameter makes everything behave as normal.

9

OK, thank you! I am currently checking if my results are better without this parameter. I will let you know, if this solves my problem. But shouldn’t a quadrature degree of 3 be OK for Langrangian polynomials of order 2 (if I recall it correctly, 9 Gauss-points are needed for the exact integration for this element type, if the integrand is a polynomial - otherwise, the quadrature will never be exact).

And can you comment on what a reasonable value for this parameter is? If FEniCS chooses the quadrature degree automatically, I have the impression, that it takes a lot of integration points which slows down the solution procedure…

10

Three sounds reasonable. In general FEniCS is quite generous with quadrature points, as it is better to over-estimate than under estimate. I would always start by using fenics automatic choice, then reduce it until you see discrepancies in the result. (Especially since you are setting this parameter globally).
You can set it locally in variational formulations through the integration measure (like: dx = Measure("dx", domain=mesh, subdomain_data=cellfunction, metadata={"quadrature_degree":2}))

11

OK, thanks a lot - for both of my problems!