Hello everyone,

I am trying to solve the following simple harmonic oscillator problem:

with this weak form:

I think I got the weak form down correctly, however, I am not getting any sensible solution, which I assume is connected to the boundary values. This is my code:

```
# %%
import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh, default_scalar_type, mesh, plot
import ufl
from ufl import dot, dx
import basix
import matplotlib.pyplot as plt
from dolfinx.fem.petsc import LinearProblem
# Define parameters
k_1 = 1.0
m_1 = 1.0
v_0 = 0.0
u_0 = 1.0
# Create mesh and define function space
domain = mesh.create_interval(MPI.COMM_WORLD, 100, (0, 5))
element = basix.ufl.element("CG", domain.basix_cell(), 1)
V = fem.functionspace(domain, basix.ufl.mixed_element([element, element]))
# Define trial and test functions
u, v = ufl.TrialFunctions(V)
u_test, v_test = ufl.TestFunctions(V)
# Define variational problem
a = (dot(u.dx(0), u_test) - dot(v, u_test)) * dx # Displacement
a += (dot(v.dx(0), v_test) + k_1/m_1 * dot(u, v_test)) * dx # Velocity
f = fem.Constant(domain, 0.0) # No driving force
L = dot(f, u_test) * dx
L += dot(f, v_test) * dx
# Apply boundary conditions
ic_facets = mesh.locate_entities_boundary(
domain, dim=0, marker=lambda x: np.isclose(x[0], 0.0)
)
dofs = fem.locate_dofs_topological(V=V.sub(0), entity_dim=0, entities=ic_facets)
bc_u = fem.dirichletbc(value=default_scalar_type(u_0), dofs=dofs, V=V.sub(0))
dofs = fem.locate_dofs_topological(V=V.sub(1), entity_dim=0, entities=ic_facets)
bc_v = fem.dirichletbc(value=default_scalar_type(v_0), dofs=dofs, V=V.sub(1))
bcs = [bc_u, bc_v]
# Solve problem
problem = LinearProblem(a, L, bcs=bcs)
solution = problem.solve()
print("Solved!")
# Split the solution to obtain the individual components
uh_u, uh_v = solution.split()
print(solution.x.array.real)
```

The code above results in an â€śinfâ€ť solution vector. Not prescribing any boundary conditions results in a zero vector solution, which makes sense if one assumes 0 initial velocity and displacement. This is why I assume that I am not prescribing the boundary conditions correctly.

Also, I am not sure if I understand the split function on the solution vector. I would assume that I would get the u and v vector respectively, however this does not seem to be the case, since the length of those arrays is the same as the solution vector (twice the number of mesh points, I would assume that it would be half of that, for u and v respectively)

Thanks in advance!