Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Re-activate commented lines in hyperelastic demo #3146

Merged
merged 4 commits into from
Apr 19, 2024
Merged

Re-activate commented lines in hyperelastic demo #3146

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
e732e51
Re-activate commented lines in hyperelastic demo
garth-wells Apr 19, 2024
e8729d2
Updates
garth-wells Apr 19, 2024
469387b
Fix constants
garth-wells Apr 19, 2024
51bdc8b
Merge branch 'main' into garth/demo-tidy
garth-wells Apr 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 6 additions & 3 deletions cpp/demo/hyperelasticity/hyperelasticity.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
from basix.ufl import element
from ufl import (
Coefficient,
Constant,
FunctionSpace,
Identity,
Mesh,
Expand All @@ -15,8 +16,10 @@
derivative,
det,
diff,
ds,
dx,
grad,
inner,
ln,
tr,
variable,
Expand All @@ -39,8 +42,8 @@

# Functions
u = Coefficient(V) # Displacement from previous iteration
# B = Coefficient(element) # Body force per unit volume
# T = Coefficient(element) # Traction force on the boundary
B = Constant(mesh, shape=(3,)) # Body force per unit volume
T = Constant(mesh, shape=(3,)) # Traction force on the boundary

# Now, we can define the kinematic quantities involved in the model:

Expand Down Expand Up @@ -72,7 +75,7 @@
psi = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2

# Total potential energy
Pi = psi * dx # - inner(B, u) * dx - inner(T, u) * ds
Pi = psi * dx - inner(B, u) * dx - inner(T, u) * ds

# First variation of Pi (directional derivative about u in the direction
# of v)
Expand Down
16 changes: 11 additions & 5 deletions cpp/demo/hyperelasticity/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,12 +156,17 @@ int main(int argc, char* argv[])
auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
functionspace_form_hyperelasticity_F_form, "u", mesh));

auto B = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
auto traction = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});

// Define solution function
auto u = std::make_shared<fem::Function<T>>(V);
auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
*form_hyperelasticity_J_form, {V, V}, {{"u", u}}, {}, {}));
auto L = std::make_shared<fem::Form<T>>(fem::create_form<T>(
*form_hyperelasticity_F_form, {V}, {{"u", u}}, {}, {}));
auto a = std::make_shared<fem::Form<T>>(
fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{"u", u}},
{{"B", B}, {"T", traction}}, {}));
auto L = std::make_shared<fem::Form<T>>(
fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{"u", u}},
{{"B", B}, {"T", traction}}, {}));

auto u_rotation = std::make_shared<fem::Function<T>>(V);
u_rotation->interpolate(
Expand Down Expand Up @@ -239,7 +244,8 @@ int main(int argc, char* argv[])
newton_solver.atol = 10 * std::numeric_limits<T>::epsilon();

la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
newton_solver.solve(_u.vec());
auto [niter, success] = newton_solver.solve(_u.vec());
std::cout << "Number of Newton iterations: " << niter << std::endl;

// Compute Cauchy stress. Construct appropriate Basix element for
// stress.
Expand Down
Loading