Skip to content

Commit

Permalink
Fix sign in von Mises demo (#16)
Browse files Browse the repository at this point in the history
* Looks reasonable.

* Fix signs.
  • Loading branch information
jhale authored Oct 31, 2024
1 parent 1127628 commit 0c0f2e7
Showing 1 changed file with 9 additions and 22 deletions.
31 changes: 9 additions & 22 deletions doc/demo/demo_plasticity_von_mises.py
Original file line number Diff line number Diff line change
Expand Up @@ -104,8 +104,8 @@
# inside the cylinder and is written as the following Neumann condition
#
# $$
# F_\text{ext}(\boldsymbol{v}) = q
# \int\limits_{\partial\Omega_\text{inner}} \boldsymbol{n} \cdot \boldsymbol{v}
# F_\text{ext}(\boldsymbol{v}) =
# \int\limits_{\partial\Omega_\text{inner}} (-q \boldsymbol{n}) \cdot \boldsymbol{v}
# \,\mathrm{d}\boldsymbol{x},
# $$
# where the vector $\boldsymbol{n}$ is the outward normal to the cylinder
Expand Down Expand Up @@ -245,7 +245,7 @@ def epsilon(v):
loading = fem.Constant(mesh, PETSc.ScalarType(0.0))

v = ufl.TestFunction(V)
F = ufl.inner(sigma, epsilon(v)) * dx - loading * ufl.inner(v, n) * ds(facet_tags_labels["inner"])
F = ufl.inner(sigma, epsilon(v)) * dx - ufl.inner(-loading*n, v) * ds(facet_tags_labels["inner"])

# Internal state
P_element = basix.ufl.quadrature_element(mesh.topology.cell_name(), degree=k_stress, value_shape=())
Expand Down Expand Up @@ -421,22 +421,6 @@ def sigma_external(derivatives):
eps = np.finfo(PETSc.ScalarType).eps
Du.x.array[:] = eps

timer = common.Timer("DOLFINx_timer")
timer.start()
evaluated_operands = evaluate_operands(F_external_operators)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
timer.stop()
pass_1 = timer.elapsed()[0]

timer.start()
evaluated_operands = evaluate_operands(F_external_operators)
_ = evaluate_external_operators(J_external_operators, evaluated_operands)
timer.stop()
pass_2 = timer.elapsed()[0]

print(f"\nNumba's JIT compilation overhead: {pass_1 - pass_2}")


# %% [markdown]
# ### Solving the problem
#
Expand All @@ -445,7 +429,7 @@ def sigma_external(derivatives):
# %%
u = fem.Function(V, name="displacement")
du = fem.Function(V, name="Newton_correction")
external_operator_problem = LinearProblem(J_replaced, F_replaced, Du, bcs=bcs)
external_operator_problem = LinearProblem(J_replaced, -F_replaced, Du, bcs=bcs)

# %%
# Defining a cell containing (Ri, 0) point, where we calculate a value of u
Expand All @@ -462,6 +446,9 @@ def sigma_external(derivatives):
loadings = q_lim * load_steps
results = np.zeros((num_increments, 2))

evaluated_operands = evaluate_operands(F_external_operators)
evaluate_external_operators(J_external_operators, evaluated_operands)

for i, loading_v in enumerate(loadings):
loading.value = loading_v
external_operator_problem.assemble_vector()
Expand All @@ -479,7 +466,7 @@ def sigma_external(derivatives):
external_operator_problem.solve(du)
du.x.scatter_forward()

Du.x.petsc_vec.axpy(-1.0, du.x.petsc_vec)
Du.x.petsc_vec.axpy(1.0, du.x.petsc_vec)
Du.x.scatter_forward()

evaluated_operands = evaluate_operands(F_external_operators)
Expand All @@ -501,7 +488,7 @@ def sigma_external(derivatives):
if MPI.COMM_WORLD.rank == 0:
print(f" it# {iteration} residual: {residual}")

u.x.petsc_vec.axpy(-1.0, Du.x.petsc_vec)
u.x.petsc_vec.axpy(1.0, Du.x.petsc_vec)
u.x.scatter_forward()

# Taking into account the history of loading
Expand Down

0 comments on commit 0c0f2e7

Please sign in to comment.