Discussion Forum

Attached is a model with one 1D element and 2 Dirichlet-type BCs on both ends of it implemented similarly to Chien Liu's blog post here: www.comsol.com/blogs/implementing-the-weak-form-in-comsol-multiphysics/

Using Linear discretisation, stationary solver, Laplace equation in the domain, this is about the simplest FEM model possible, so it should be very easy to interpret.

As described in Chien Liu's blog post that follows the one linked above in that series
www.comsol.com/blogs/implementing-the-weak-form-with-a-comsol-app
we can add an Assemble node under the solver node, and view the System Matrix under Derived Values.

In our case we have 4 DOFs: two corresponding to the dependent variable's value at each end of the element, and two for the lambdas used to implement the Dirichlet BCs. The stiffness matrix looks about right:
0 1 0 0
1 1 0 -1
0 0 0 1
0 -1 1 1
But the load vector is all zeros:
0
0
0
0
(Sometimes when testing with some trivial changes, the middle entries show up as small values around the machine epsilon ~1e-16.)

The BCs are 2.0 on the left end and 7.0 on the right, and the solution shows up correctly in the plot.

The eliminated stiffness matrix and load vector are the same (up to a scaling constant), because we are using weak constraints (nothing is eliminated, and the Lagrange multipliers are explicitly solved for). If we implement the Dirichlet BCs the usual way instead, and check the "Use weak constraints" checkbox, then exactly the same thing happens as above. Also, we can see under Equation View, that we have the same equations as when we use the "Weak Contribution" method to implement the Dirichlet BCs.

If we use LiveLink and get the Kc, Lc matrices, then, of course, the solution Kc\Lc is just the trivial all-zeros solution (as Lc is a zero vector, and so is L), which is incorrect.

There is a comment under Chien Liu's post (second link above) from Jing Wang July 6, 2015, reporting the same issue, which has been referred to Comsol support, but with no public follow-up. My guess is that, since the blog post refers to Comsol 5.0, that it used to work as expected under 5.0, but has changed since then.

The interesting thing is that I found a way to get the expected load vector, using this strange procedure: in "Stationary Solver 1" under "Linearity" choose "Linear perturbation", with "Values of linearization point" "Prescribed by" "Solution", "Solution": "Solution 1 (sol1)". If we solve and assemble now, we get a flat zero solution, but the load vector is now evaluated to the expected values.

However, under no combination of these Linearity settings could I get both the expected solution in the plots and the expected load vector from System Matrix evaluation.

My interpretation is that either this is a bug, or Comsol since 5.1 has been using another layer of linearisation point picking, that we can't access. Either way, we can't export the matrices into Matlab and solve the system there with higher precision using the procedure dictated by common sense.