-
Notifications
You must be signed in to change notification settings - Fork 3
/
dummy.py
58 lines (41 loc) · 1.93 KB
/
dummy.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import numpy as np
import time
import precice
# preCICE setup
participant_name = "Dummy"
config_file_name = "precice-config.xml"
solver_process_index = 0
solver_process_size = 1
interface = precice.Interface(participant_name, config_file_name, solver_process_index, solver_process_size)
mesh_name = "Dummy-Mesh"
mesh_id = interface.get_mesh_id(mesh_name)
lmbda_id = interface.get_data_id("Lmbda", mesh_id)
mu_id = interface.get_data_id("Mu", mesh_id)
gc_id = interface.get_data_id("Gc", mesh_id)
# define coupling mesh
L = 1e-3 # domain size in m
N = 100 # number of cells in each direction
axis = np.linspace(-L/2,L/2,N+1) # coordinates along one axis
positions = np.transpose([np.tile(axis, len(axis)), np.repeat(axis, len(axis))]) # mesh coordinates
vertex_ids = interface.set_mesh_vertices(mesh_id, positions)
precice_dt = interface.initialize() # pseudo timestep size handled by preCICE
step = 0
while interface.is_coupling_ongoing():
print("Generating data")
time.sleep(1.0) # for better readability of shell output
lmbda = np.ones((N+1)*(N+1)) * 121153.8e6 # First Lamé parameter in Pa, constant everywhere
mu = np.ones((N+1)*(N+1)) * 80769.2e6 # Second Lamé parameter in Pa, constant everywhere
gc = np.ones((N+1)*(N+1)) * 2.7e3 # Fracture toughness in N/m, constant everywhere
# modify fracture toughness at some arbitrary location (here vertical line, verically centered, right quarter of domain)
for i in range(int(N*3/4), N+1):
gc[(N+1)*int(N/2) + i] *= np.exp(-step) # fracture toughness decreases exponentially
gc[(N+1)*(int(N/2)+1)+ + i] *= np.exp(-step)
print(np.exp(-step))
# write data to preCICE
interface.write_block_scalar_data(lmbda_id, vertex_ids, lmbda)
interface.write_block_scalar_data(mu_id, vertex_ids, mu)
interface.write_block_scalar_data(gc_id, vertex_ids, gc)
# do the coupling
precice_dt = interface.advance(precice_dt)
step = step + 1
interface.finalize()