Plotting graphs

General
1

hi all,

im trying to plot a graph of stress against deformation when applying force on a rod.
what im trying to do it to create a loop that is applying on each loop more force on the rod until i reach my max force. i want in the graph to see the avg deformation at each time point, how can i do it? is there a specific method to achieve that? let me know if you need any more info for that

(since its non-homogenous hyperplastic material i dont want to plot it on one point only, i want the avg deformation on all points or specific ones)

2

Dear @Mirialex,

You have not specified what what version of FEniCS your are using (legacy: what version, DOLFINx: what version).

You have also not specified if you want the average magnitude of deformation, or the average in each coordinate direction (x,y,z).

Please answer these questions to increase the likelihood of getting a useful reply

3

hi,
this is the version i have 2019.2.0.dev0
and also i want average deformation for (x,0.1,0.1) for example at each time point . the average deformation for all Xs on the line of y=z=0.1 for each time step.
or perhaps the average magnitude of deformation at each time step in the loop

4

Does your mesh align with this line?

5

yes sure, this is my mesh
length, width, height = 1.0, 0.2, 0.2 # Dimensions of the rod
num_elements = 40 # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 5, 5)

6

You will have issues with integrating over a line, which is an edge of the mesh (and the line is not part of the mesh).

I will illustrate how to do the average over a plane (y=0.1) which aligns with the grid (note that I changed number of elements from 5 to 4 to make this work).

from dolfin import *


class IntegralBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.1)


length, width, height = 1.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 40  # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 4, 4)

ff = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
IntegralBoundary().mark(ff, 1)


dS_plane = Measure("dS", domain=mesh, subdomain_data=ff, subdomain_id=1)

V = VectorFunctionSpace(mesh, "Lagrange", 1)
u = Function(V)
u.interpolate(Expression(("x[0]", "x[1]", "x[2]"), degree=1))

vol_plane = assemble(Constant(1) * dS_plane)
int_u = assemble(sqrt(dot(u, u)) * dS_plane)
print(int_u / vol_plane)

If you want to do the average over a line not part of your mesh, you should get a quadrature rule for the interval (of your choice), and use

from dolfin import *
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np

shape = "interval"
deg = 10
scheme = "default"

points, weights = cquad(shape, deg, scheme)


x_start = -1.0
length, width, height = 2.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 40  # Number of elements along each dimension
mesh = BoxMesh(Point(x_start, 0, 0), Point(length, width, height), num_elements, 4, 4)


V = VectorFunctionSpace(mesh, "Lagrange", 1)
u = Function(V)
u.interpolate(Expression(("x[0]", "x[1]", "x[2]"), degree=1))


absJ = abs(length)
physical_points = np.zeros((len(points), 3))
physical_points[:, 0] = np.asarray([x_start + p * absJ for p in points]).reshape(-1)
physical_points[:, 1] = 0.1
physical_points[:, 2] = 0.2

values = [np.sqrt(np.dot(u(point), u(point))) for point in physical_points]
compute_integral = np.dot(values, weights) * absJ

print(compute_integral)
# Exact value is `int_-1^1 sqrt(x**2 + 0.1*0.1 + 0.2*0.2) dx = 1.13486`
# https://www.wolframalpha.com/input?i=int_-1%5E1+sqrt%28x**2+%2B+0.1*0.1+%2B+0.2*0.2%29+dx

yielding

1.1382180980566456

with 10 as quadrature degree. You can reduce the precision to observe a reduction in the accuracy.

7

ok i will check on that thanks, what do you mean its not part of the mesh? and also why its not part of the mesh?
not all the points that are between 0-0.2 are part of the mesh?

8

also, i dont understand how from here i create a graph plot of the stress and deformation?

9

They are inside thw mesh, but the line is not aligned with any edges of the mesh.

You asked about how to plot the average stress along a line in a loop where you apply more and more force.

I have shown how to get the average u along a ling for a given load. Put this inside a for-loop, and store the results in an array, and use matplotlib for plotting the result.

10

hi, thanks for the reply
couple of questions questions:

  1. why do you start the mesh from -1?
    "x_start = -1.0
    length, width, height = 2.0, 0.2, 0.2 # Dimensions of the rod
    num_elements = 40 # Number of elements along each dimension
    mesh = BoxMesh(Point(x_start, 0, 0), Point(length, width, height), num_elements, 4, 4)
    "
  2. im getting " *** Warning: Found no facets matching domain for boundary condition."

those are my boundary conditions

x_start = -1.0
length, width, height = 2.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 40  # Number of elements along each dimension
mesh = BoxMesh(Point(x_start, 0, 0), Point(length, width, height), num_elements, 4, 4)

# mesh = UnitCubeMesh(24, 16, 16)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

#creating subtomains and boundaty conditions

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 1 + DOLFIN_EPS


class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > 1 - DOLFIN_EPS


subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)

boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

left = AutoSubDomain(lambda x: near(x[0], 0))
right = AutoSubDomain(lambda x: near(x[0], 2))
up = AutoSubDomain(lambda x: near(x[2], 0.2))
down = AutoSubDomain(lambda x: near(x[2], 0))
front = AutoSubDomain(lambda x: near(x[1], 0))
back = AutoSubDomain(lambda x: near(x[1], 0.2))

left.mark(boundary_parts, 1)
right.mark(boundary_parts, 2)
up.mark(boundary_parts, 5)
down.mark(boundary_parts, 6)
front.mark(boundary_parts, 3)
back.mark(boundary_parts, 4)

bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 1)
bc3 = DirichletBC(V.sub(1), Constant(0), boundary_parts, 3)
bc5 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 5)
bc6 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 6)
tDirBC = Expression(('1.0*time_'), time_=0.0, degree=0)
tDirBC4 = Expression(('0.1*time_'), time_=0.0, degree=0)
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 2)
bc4 = DirichletBC(V.sub(1), tDirBC4, boundary_parts, 4)
#bc2 = DirichletBC(V.sub(0), Constant(1), boundary_parts, 2)
#bc4 = DirichletBC(V.sub(1), Constant(0.2), boundary_parts, 4)
# Condition for one side

bcs = [bc1, bc2, bc3, bc4, bc5, bc6]

**ps, managed to over come the error when i changed left = AutoSubDomain(lambda x: near(x[0], 0)) to left = AutoSubDomain(lambda x: near(x[0], -1))

11

also, if i use the second option compute_integral = np.dot(values, weights) * absJ

do i need to use it as append in my for loop ?
cause if i use it with append i get an array of arrays
for exp : [array([-4.04268259e-11, -8.51276143e-11, -1.10411973e-10, -1.10411973e-10,
-8.51276143e-11, -4.04268259e-11]), array([7.42936737e-11, 1.56441745e-10, 2.02907621e-10, 2.02907621e-10,
1.56441745e-10, 7.42936737e-11]), array([2.96305786e-10, 6.23937303e-10, 8.09257360e-10, 8.09257360e-10,
6.23937303e-10, 2.96305786e-10]), array([-4.03758192e-10, -8.50202083e-10, -1.10272665e-09, -1.10272665e-09,
-8.50202083e-10, -4.03758192e-10]), array([1.40439230e-11, 2.95725829e-11, 3.83561461e-11, 3.83561461e-11,
2.95725829e-11, 1.40439230e-11]), array([-6.01814288e-11, -1.26725296e-10, -1.64364877e-10, -1.64364877e-10,
-1.26725296e-10, -6.01814288e-11]), array([2.14710184e-10, 4.52119732e-10, 5.86407030e-10, 5.86407030e-10,
4.52119732e-10, 2.14710184e-10]), array([-1.66728297e-10, -3.51083268e-10, -4.55361008e-10, -4.55361008e-10,
-3.51083268e-10, -1.66728297e-10]), array([1.28421937e-11, 2.70420763e-11, 3.50740358e-11, 3.50740358e-11,
2.70420763e-11, 1.28421937e-11]), array([-9.03286021e-11, -1.90206830e-10, -2.46701514e-10, -2.46701514e-10,
-1.90206830e-10, -9.03286021e-11]), array([-5.05214331e-11, -1.06384040e-10, -1.37981921e-10, -1.37981921e-10,
-1.06384040e-10, -5.05214331e-11]), array([-2.88591664e-11, -6.07693514e-11, -7.88188888e-11, -7.88188888e-11,
-6.07693514e-11, -2.88591664e-11]), array([-1.34599776e-11, -2.83429568e-11, -3.67613001e-11, -3.67613001e-11,
-2.83429568e-11, -1.34599776e-11]), array([-5.67405218e-11, -1.19479705e-10, -1.54967223e-10, -1.54967223e-10,
-1.19479705e-10, -5.67405218e-11])]
[array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449]), array([0.17132449, 0.36076157, 0.46791393, 0.46791393, 0.36076157,
0.17132449])]
one for the deformation and one for the stress

12

This was just to illustrate that it works for a mesh with start at non-zero. Feel free to change it back (that might resolve your boundary condition issues as well).

13

Compute integral should return a single scalar value. Please make a minimal reproducible example that yields the behavior above.

14

this is the loop im currently using withot the append that was causing what i have sent.
but i worry that it just takes the last result of the loop and than plot that last result.

stress_values_all = []
deformation_gradient_values_all = []
stress_values_all1 = []
deformation_gradient_values_all2 = []
stress_values_all2 = []
deformation_gradient_values_all3 = []

stress123 = Function(S, name='stressi')
dt = 0.1
t, T = 0.0, 100 * dt
TT = Traction()

#creating functions in order to see the stress and deformation over time
V_tensor = TensorFunctionSpace(mesh, "P", 1)
defGrad = Function(V_tensor, name='F')
# F_projected = project(F, V_tensor)
# F_xx = project(F_projected[0, 0], S)
shape = "interval"
deg = 10
scheme = "default"

points, weights = cquad(shape, deg, scheme)

#creating a file for paraview
stress123.rename("SRESS1", "")
materials_function.rename("materials", "")
defGrad.rename("Deformation_DRAD", "")
u.rename("Displacement Vector", "")
test_file = fe.XDMFFile("test.xdmf")
test_file.parameters["flush_output"] = True
test_file.parameters["functions_share_mesh"] = True

absJ = abs(length)
physical_points = np.zeros((len(points), 3))
physical_points[:, 0] = np.asarray([x_start + p * absJ for p in points]).reshape(-1)
physical_points[:, 1] = 0.1
physical_points[:, 2] = 0.1
u.rename("u", "displacement")
#creating a loop with time dependent in order to apply force over time
while t <= T:
    print('time: ', t)

    # Increase traction
    TT.t = t
    tDirBC.time_ = t

    solve(L == 0, u, bcs, solver_parameters={"newton_solver": {"linear_solver": "lu",
                                                               "relative_tolerance": 1e-6,
                                                               # "preconditioner":"ilu",
                                                               "convergence_criterion": "incremental",
                                                               }})

    # get stretch at a point for plotting
    #trying to create avarage stress and deformation

    for point in physical_points:
        # is to return the coordinates of all degrees of freedom (DOFs) associated with the finite element function space "S."

        try:
            z = 0
            # is to return the coordinates of all degrees of freedom (DOFs) associated with the finite element function space "S."
            DF = I + grad(u)  # Compute the deformation gradient
            defGrad.assign(project(DF, V_tensor))  # Project the deformation gradient
            deformation_gradient_values_all=defGrad(point)[0]  # Store the value at the current point

            sigma_xx = sigma[0, 0]  # Extravgacting the xx component of the stress tensor
            stress_xx_projected = project(sigma_xx, S)  # Projecting this scalar component
            stress123.assign(stress_xx_projected)  # Assigning the projected stress
            stress_values_all=stress123(point)  # Store the stress value at the current point

        except Exception as e:
            print(f"Could not evaluate at {point}: {e}")

    stress_values_all1=(np.dot(stress_values_all, weights) * absJ)
    deformation_gradient_values_all2=(np.dot(deformation_gradient_values_all, weights) * absJ)
    print(stress_values_all1)
    print(deformation_gradient_values_all2)
    print(stress_values_all)
    print(deformation_gradient_values_all)



    # save xdmf file
    test_file.write(u, t)
    test_file.write(materials_function, t)
    test_file.write(stress123, t)
    test_file.write(defGrad, t)

    # time increment
    t += float(dt)

deformation_gradient_values_all3 = np.array(deformation_gradient_values_all2)
stress_values_all2 = np.array(stress_values_all1)
15

You are not appending to a list, so thus you cannot expect it to store all values. Please make small modifications to my original minimal example, to understand how the code would work. I do not have the bandwidth to consider non-reprodicble code (i.e. code I cannot copy paste into an editor and run to get an error message). Also you should first understand how to compute the average displacement magnitude in a loop before moving to stress computations

16

this is my printing after trying with append when trying to use your exp for the deformation

"
[array([-2.99186406e-12, -6.04402457e-12, -7.18381208e-12, -6.04402457e-12,
-2.99186406e-12]), array([20036.33126886, 40476.46426713, 48109.55173048, 40476.46426713,
20036.33126886]), array([39991.46538114, 80788.89781613, 96024.13968479, 80788.89781613,
39991.46538114]), array([ 59893.44315511, 120993.94740349, 143811.09311492, 120993.94740349,
59893.44315511]), array([ 79758.60863172, 161124.63050036, 191509.65595594, 161124.63050036,
79758.60863172]), array([ 99597.13805365, 201201.50467878, 239144.26254503, 201201.50467878,
99597.13805365]), array([119415.70405786, 241238.04968948, 286730.93465627, 241238.04968948,
119415.70405786])]
[array([6.12416715e-20, 1.23717575e-19, 1.47048345e-19, 1.23717575e-19,
6.12416715e-20]), array([0.0108448 , 0.02190816, 0.02603962, 0.02190816, 0.0108448 ]), array([0.02168924, 0.04381558, 0.05207837, 0.04381558, 0.02168924]), array([0.03253345, 0.06572257, 0.07811659, 0.06572257, 0.03253345]), array([0.04337754, 0.08762929, 0.1041545 , 0.08762929, 0.04337754]), array([0.05422155, 0.10953585, 0.13019222, 0.10953585, 0.05422155]), array([0.06506551, 0.13144232, 0.15622983, 0.13144232, 0.06506551])]
"

this is the loop i created


# List to store stress and deformation gradient values for all points
stress_values_all = []
deformation_gradient_values_all = []
stress_values_all1 = []
deformation_gradient_values_all2 = []
stress_values_all2 = []
deformation_gradient_values_all3 = []

stress123 = Function(S, name='stressi')
dt = 0.1
t, T = 0.0, 50 * dt
TT = Traction()

#creating functions in order to see the stress and deformation over time
V_tensor = TensorFunctionSpace(mesh, "P", 1)
defGrad = Function(V_tensor, name='F')
stress123 = Function(V_tensor, name='stressi')
# F_projected = project(F, V_tensor)
# F_xx = project(F_projected[0, 0], S)
shape = "interval"
deg = 8
scheme = "default"

points, weights = cquad(shape, deg, scheme)

#creating a file for paraview
stress123.rename("SRESS1", "")
materials_function.rename("materials", "")
defGrad.rename("Deformation_DRAD", "")
u.rename("Displacement Vector", "")
test_file = fe.XDMFFile("test.xdmf")
test_file.parameters["flush_output"] = True
test_file.parameters["functions_share_mesh"] = True

absJ = abs(length)
physical_points = np.zeros((len(points), 3))
physical_points[:, 0] = np.asarray([p * absJ for p in points]).reshape(-1)
physical_points[:, 1] = 0.2
physical_points[:, 2] = 0.2
u.rename("u", "displacement")
#creating a loop with time dependent in order to apply force over time
while t <= T:
    print('time: ', t)

    # Increase traction
    TT.t = t
    tDirBC.time_ = t

    solve(L == 0, u, bcs, solver_parameters={"newton_solver": {"linear_solver": "lu",
                                                               "relative_tolerance": 1e-6,
                                                               # "preconditioner":"ilu",
                                                               "convergence_criterion": "incremental",
                                                               }})

    # get stretch at a point for plotting
    #trying to create avarage stress and deformation

    for point in physical_points:
        # is to return the coordinates of all degrees of freedom (DOFs) associated with the finite element function space "S."

        try:
            z = 0
            # is to return the coordinates of all degrees of freedom (DOFs) associated with the finite element function space "S."
            DF = I + grad(u)  # Compute the deformation gradient
            defGrad.assign(project(DF, V_tensor))  # Project the deformation gradient
            #deformation_gradient_values_all=defGrad(point)[0]  # Store the value at the current point
            deformation_gradient_values_all = [np.sqrt(np.dot(u(point), u(point)))][0]

            sigma_xx = sigma# Extravgacting the xx component of the stress tensor
            stress_xx_projected = project(sigma_xx, V_tensor)  # Projecting this scalar component
            stress123.assign(stress_xx_projected)  # Assigning the projected stress
            stress_values_all=stress123(point)[0]  # Store the stress value at the current point

        except Exception as e:
            print(f"Could not evaluate at {point}: {e}")

    stress_values_all1.append(np.dot(stress_values_all, weights) * absJ)
    deformation_gradient_values_all2.append(np.dot(deformation_gradient_values_all, weights) * absJ)
    print(stress_values_all1)
    print(deformation_gradient_values_all2)
    print(stress_values_all)
    print(deformation_gradient_values_all)

    # save xdmf file
    test_file.write(u, t)
    test_file.write(materials_function, t)
    test_file.write(stress123, t)
    test_file.write(defGrad, t)

    # time increment
    t += float(dt)

deformation_gradient_values_all3 = np.array(deformation_gradient_values_all2)
stress_values_all2 = np.array(stress_values_all1)

any idea now why its creating an array of arrays?

**ps
managed to overcome the issue with the array of arrays after setting append in here stress_values_all=stress123(point)[0]

17

Why is the dot product outside the loop?

18

i have entered it to the loop

another question, if i would want to get the sigma on y and z, whould it be right to change it to 1 or to accordantly in the loop? " stress_values_all.append(stress123(point)[2])"

            deformation_gradient_values_all.append(defGrad(point)[2])# Store the value at the current point

            stress_values_all.append(stress123(point)[2])
            stress_values_all1.append(np.dot(stress_values_all, weights) * absJ)
            deformation_gradient_values_all2.append(np.dot(deformation_gradient_values_all, weights) * absJ)

this is how stress is defined

V_tensor = TensorFunctionSpace(mesh, "P", 1)
defGrad = Function(V_tensor, name='F')
stress123 = Function(V_tensor, name='stressi')

   sigma = (1 / J) * dev_stress + vol_stress
    u.rename("u", "displacement")

    DF = I + grad(u)  # Compute the deformation gradient
    defGrad.assign(project(DF, V_tensor))  # Project the deformation gradient
    sigma_xx = sigma # Extravgacting the xx component of the stress tensor
    stress_xx_projected = project(sigma_xx, V_tensor)  # Projecting this scalar component
    stress123.assign(stress_xx_projected)
19

Seems to me that you decided to ask the question from the last post in a new topic, at Extracting sigmayy and sigmazz

In that case, I’ll close this one. Get in touch with me if I misunderstood.

1 Like