Skip to content

Commit

Permalink
WIP: test for example 3 and 1
Browse files Browse the repository at this point in the history
  • Loading branch information
@OmarDuran
OmarDuran committed Sep 18, 2024
1 parent cbbaf53 commit 8b31524
Show file tree
Hide file tree
Showing 3 changed files with 119 additions and 76 deletions.
19 changes: 11 additions & 8 deletions examples/cosserat_elasticity/compute_results.sh
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,16 @@
#!/bin/bash


echo " "
echo "Example 3:: started"
python -OO mfe_cosserat_elasticity_example_3.py
mkdir output_example_3
mv wc_rt* output_example_3/
mv wc_bdm* output_example_3/
mv geometric_* output_example_3/
echo "Example 3:: completed"


echo " "
echo "Example 1:: started"
python -OO mfe_cosserat_elasticity_example_1.py
Expand All @@ -25,14 +35,7 @@ echo "Example 1:: completed"
# echo "Example 2:: completed"


# echo " "
# echo "Example 3:: started"
# python -OO mfe_cosserat_elasticity_example_3.py
# mkdir output_example_3
# mv wc_rt* output_example_3/
# mv wc_bdm* output_example_3/
# mv geometric_* output_example_3/
# echo "Example 3:: completed"




Expand Down
93 changes: 52 additions & 41 deletions examples/cosserat_elasticity/mfe_cosserat_elasticity_example_2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,19 +78,26 @@ def create_product_space(method, gmesh):
"t": t_disc_Q,
}

physical_tags = {
"s": [1],
"m": [1],
"u": [1],
"t": [1],
}

s_field_bc_physical_tags = [2, 3, 4, 5]
m_field_bc_physical_tags = [2, 3, 4, 5]
if gmesh.dimension == 3:
s_field_bc_physical_tags = [2, 3, 4, 5, 6, 7]
m_field_bc_physical_tags = [2, 3, 4, 5, 6, 7]
discrete_spaces_bc_physical_tags = {
b_physical_tags = {
"s": s_field_bc_physical_tags,
"m": m_field_bc_physical_tags,
}

space = ProductSpace(discrete_spaces_data)
space.make_subspaces_discontinuous(discrete_spaces_disc)
space.build_structures(discrete_spaces_bc_physical_tags)
space.build_structures(physical_tags, b_physical_tags)
return space


Expand Down Expand Up @@ -119,39 +126,16 @@ def four_field_approximation(material_data, method, gmesh, symmetric_solver_q=Tr
P.setType("sbaij")

# Material data
m_lambda = material_data["lambda"]
m_mu = material_data["mu"]
m_kappa = material_data["kappa"]
m_gamma = material_data["gamma"]

# exact solution
u_exact = lce.displacement(m_lambda, m_mu, m_kappa, m_gamma, dim)
t_exact = lce.rotation(m_lambda, m_mu, m_kappa, m_gamma, dim)
s_exact = lce.stress(m_lambda, m_mu, m_kappa, m_gamma, dim)
m_exact = lce.couple_stress(m_lambda, m_mu, m_kappa, m_gamma, dim)
div_s_exact = lce.stress_divergence(m_lambda, m_mu, m_kappa, m_gamma, dim)
div_m_exact = lce.couple_stress_divergence(m_lambda, m_mu, m_kappa, m_gamma, dim)
f_rhs = lce.rhs(m_lambda, m_mu, m_kappa, m_gamma, dim)

def f_lambda(x, y, z):
return m_lambda

def f_mu(x, y, z):
return m_mu

def f_kappa(x, y, z):
return m_kappa

def f_gamma(x, y, z):
return m_gamma
u_exact = lce.displacement(material_data, dim)
t_exact = lce.rotation(material_data, dim)
s_exact = lce.stress(material_data, dim)
m_exact = lce.couple_stress(material_data, dim)
div_s_exact = lce.stress_divergence(material_data, dim)
div_m_exact = lce.couple_stress_divergence(material_data, dim)
f_rhs = lce.rhs(material_data, dim)

m_functions = {
"rhs": f_rhs,
"lambda": f_lambda,
"mu": f_mu,
"kappa": f_kappa,
"gamma": f_gamma,
}
m_functions = lce.get_material_functions(material_data, dim)
m_functions["rhs"] = f_rhs

exact_functions = {
"s": s_exact,
Expand Down Expand Up @@ -736,29 +720,56 @@ def method_definition(k_order):

methods = [method_1, method_2, method_3, method_4]
method_names = ["sc_rt", "sc_bdm", "wc_rt", "wc_bdm"]

methods = [method_3]
method_names = ["wc_rt"]
return zip(method_names, methods)


def material_data_definition():
# Material data for example 2
case_0 = {"lambda": 1.0, "mu": 1.0, "kappa": 1.0, "gamma": 1.0}
case_1 = {"lambda": 1.0e2, "mu": 1.0, "kappa": 1.0, "gamma": 1.0}
case_2 = {"lambda": 1.0e4, "mu": 1.0, "kappa": 1.0, "gamma": 1.0}
case_0 = {
"lambda_s": 1.0,
"mu_s": 1.0,
"kappa_s": 0.1,
"lambda_o": 1.0,
"mu_o": 1.0,
"kappa_o": 0.1,
"l": 1.0,
}
case_1 = {
"lambda_s": 1.0e2,
"mu_s": 1.0,
"kappa_s": 0.1,
"lambda_o": 1.0,
"mu_o": 1.0,
"kappa_o": 0.1,
"l": 1.0,
}
case_2 = {
"lambda_s": 1.0e4,
"mu_s": 1.0,
"kappa_s": 0.1,
"lambda_o": 1.0,
"mu_o": 1.0,
"kappa_o": 0.1,
"l": 1.0,
}
cases = [case_0, case_1, case_2]

cases = [case_2]
return cases


def main():
approximation_q = True
postprocessing_q = True
refinements = {0: 4, 1: 4}
refinements = {0: 2, 1: 4}
case_data = material_data_definition()
for k in [0, 1]:
for k in [0]:
methods = method_definition(k)
for i, method in enumerate(methods):
for material_data in case_data:
for d in [3]:
for d in [2]:
configuration = {
"k_order": k,
"dimension": d,
Expand Down
83 changes: 56 additions & 27 deletions examples/cosserat_elasticity/strong_solution_cosserat_elasticity_example_1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -397,72 +397,101 @@ def couple_stress_scaled(material_data, dim: int = 2):
]
)
else:
return lambda x, y, z: m_l * np.array(
return lambda x, y, z: np.array(
[
[
-(
(-1 + 2 * x)
* (m_lambda_o + 2 * m_mu_o)
* np.sin(np.pi * y)
* np.sin(np.pi * z)
)
+ m_lambda_o
* np.sin(np.pi * x)
* (
(1 - 2 * z) * np.sin(np.pi * y)
+ (1 - 2 * y) * np.sin(np.pi * z)
m_l
* (
(-1 + 2 * x)
* (m_lambda_o + 2 * m_mu_o)
* np.sin(np.pi * y)
* np.sin(np.pi * z)
+ m_lambda_o
* np.sin(np.pi * x)
* (
(-1 + 2 * z) * np.sin(np.pi * y)
+ (-1 + 2 * y) * np.sin(np.pi * z)
)
)
),
np.pi
m_l
* np.pi
* (
(-1 + y) * y * (m_kappa_o - m_mu_o) * np.cos(np.pi * x)
- (-1 + x) * x * (m_kappa_o + m_mu_o) * np.cos(np.pi * y)
)
* np.sin(np.pi * z),
np.pi
m_l
* np.pi
* (
(-1 + z) * z * (m_kappa_o - m_mu_o) * np.cos(np.pi * x)
- (-1 + x) * x * (m_kappa_o + m_mu_o) * np.cos(np.pi * z)
)
* np.sin(np.pi * y),
],
[
np.pi
m_l
* np.pi
* (
-((-1 + y) * y * (m_kappa_o + m_mu_o) * np.cos(np.pi * x))
+ (-1 + x) * x * (m_kappa_o - m_mu_o) * np.cos(np.pi * y)
)
* np.sin(np.pi * z),
(1 - 2 * x) * m_lambda_o * np.sin(np.pi * y) * np.sin(np.pi * z)
+ np.sin(np.pi * x)
* (
(m_lambda_o - 2 * z * m_lambda_o) * np.sin(np.pi * y)
- (-1 + 2 * y) * (m_lambda_o + 2 * m_mu_o) * np.sin(np.pi * z)
-(
m_l
* (
(-1 + 2 * x)
* m_lambda_o
* np.sin(np.pi * y)
* np.sin(np.pi * z)
+ np.sin(np.pi * x)
* (
(-1 + 2 * z) * m_lambda_o * np.sin(np.pi * y)
+ (-1 + 2 * y)
* (m_lambda_o + 2 * m_mu_o)
* np.sin(np.pi * z)
)
)
),
np.pi
m_l
* np.pi
* (
(-1 + z) * z * (m_kappa_o - m_mu_o) * np.cos(np.pi * y)
- (-1 + y) * y * (m_kappa_o + m_mu_o) * np.cos(np.pi * z)
)
* np.sin(np.pi * x),
],
[
np.pi
m_l
* np.pi
* (
-((-1 + z) * z * (m_kappa_o + m_mu_o) * np.cos(np.pi * x))
+ (-1 + x) * x * (m_kappa_o - m_mu_o) * np.cos(np.pi * z)
)
* np.sin(np.pi * y),
np.pi
m_l
* np.pi
* (
-((-1 + z) * z * (m_kappa_o + m_mu_o) * np.cos(np.pi * y))
+ (-1 + y) * y * (m_kappa_o - m_mu_o) * np.cos(np.pi * z)
)
* np.sin(np.pi * x),
(1 - 2 * x) * m_lambda_o * np.sin(np.pi * y) * np.sin(np.pi * z)
+ np.sin(np.pi * x)
* (
-((-1 + 2 * z) * (m_lambda_o + 2 * m_mu_o) * np.sin(np.pi * y))
+ (1 - 2 * y) * m_lambda_o * np.sin(np.pi * z)
-(
m_l
* (
(-1 + 2 * x)
* m_lambda_o
* np.sin(np.pi * y)
* np.sin(np.pi * z)
+ np.sin(np.pi * x)
* (
(-1 + 2 * z)
* (m_lambda_o + 2 * m_mu_o)
* np.sin(np.pi * y)
+ (-1 + 2 * y) * m_lambda_o * np.sin(np.pi * z)
)
)
),
],
]
Expand Down

0 comments on commit 8b31524

Please sign in to comment.