|
46 | 46 | domain = mesh.create_box(MPI.COMM_WORLD, extent, [Nx, Ny, Nz]) |
47 | 47 |
|
48 | 48 | # create the function space |
49 | | -V = fem.FunctionSpace(domain, ("Lagrange", 1, (domain.geometry.dim,))) |
| 49 | +V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,))) |
50 | 50 |
|
51 | 51 | # define some tags |
52 | 52 | tag_left, tag_top, tag_right, tag_bottom, tag_back, tag_front = 1, 2, 3, 4, 5, 6 |
|
152 | 152 | kwargsTScheme = dict(scheme='g-a-newmark', alpha_m=alpha_m, alpha_f=alpha_f) |
153 | 153 |
|
154 | 154 | # Time integration: define a TimeStepper instance |
155 | | -tStepper = TimeStepper.build(V, pde.m, pde.c, pde.k, pde.L, dt, bcs=bcs, **kwargsTScheme) |
| 155 | +tStepper = TimeStepper.build(V, |
| 156 | + pde.M_fn, pde.C_fn, pde.K_fn, pde.b_fn, dt, bcs=bcs, |
| 157 | + **kwargsTScheme) |
156 | 158 |
|
157 | 159 | # Set the initial values |
158 | 160 | tStepper.set_initial_condition(u0=[0, 0, 0], v0=[0, 0, 0], t0=0) |
|
184 | 186 | u_n = tStepper.timescheme.u # The displacement at time t_n |
185 | 187 | v_n = tStepper.timescheme.v # The velocity at time t_n |
186 | 188 | comm = V.mesh.comm |
187 | | -Energy_elastic = lambda *a: comm.allreduce(fem.assemble_scalar(fem.form(1/2 * pde.k(u_n, u_n))), op=MPI.SUM) |
188 | | -Energy_kinetic = lambda *a: comm.allreduce(fem.assemble_scalar(fem.form(1/2 * pde.m(v_n, v_n))), op=MPI.SUM) |
189 | | -Energy_damping = lambda *a: dt * comm.allreduce(fem.assemble_scalar(fem.form(pde.c(v_n, v_n))), op=MPI.SUM) |
| 189 | + |
| 190 | + |
| 191 | +def Energy_elastic(): |
| 192 | + return comm.allreduce(fem.assemble_scalar(fem.form(1/2 * pde.K_fn(u_n, u_n))), op=MPI.SUM) |
| 193 | + |
| 194 | + |
| 195 | +def Energy_kinetic(): |
| 196 | + return comm.allreduce(fem.assemble_scalar(fem.form(1/2 * pde.M_fn(v_n, v_n))), op=MPI.SUM) |
| 197 | + |
| 198 | + |
| 199 | +def Energy_damping(): |
| 200 | + return dt * comm.allreduce(fem.assemble_scalar(fem.form(pde.C_fn(v_n, v_n))), op=MPI.SUM) |
| 201 | + |
190 | 202 |
|
191 | 203 | # -> live plotting parameters |
192 | 204 | clim = 0.4 * L_ * B_ * H_ / (E * B_ * H_**3 / 12) * np.amax(F_0) * np.array([0, 1]) |
|
0 commit comments