DG and Newtons method linearization

Dear All,

First of all, congratulations to the developers of Netgen. I think it is a very good solver and what is more, the documentation is also properly developed which is not that usual.

I’m working on a solver for a system of 3 nonlinear equations. Following the tutorials and documentation, I wrote a Newton solver that solves the problem using a 1D mesh. As some of the parameters in my equations are discontinuous I got into troubles and wanted to implement a DG method to include the discontinuity of the parameters. I followed the DDG tutorial and added the jumps on the boundaries of the elements. In my case that would be points as I am working with 1D mesh.

Unfortunately, during the linearization of the form, I got an error

RuntimeError: CalcLinearization for element-wise skeleton-VOL not implementedin AssembleLinearization

Is there anything I can do to overcome this problem or is it just something that is not implemented and won’t change?


could you please attach your python file.
Sounds like you used a combination of arguments in the SymbolicBFI which is not implemented.

Best regards


Here is the interesting part of the code. I cut it out of a bigger code to show where the problem is.

Best wishes,

Attachment: PNP.py

Hi B,

thank you for the flowers …

it is as the error tells us, linearization for DG is not implemented:
RuntimeError: CalcLinearization for element-wise skeleton-VOL not implementedin AssembleLinearization

you can replace your DG formulation by an hybrid-DG (HDG), then it shall work.

It is not a principal problem with non-linear DG, but we mainly use HDG (rather than DG), and so nobody found time to implement that.

Best, Joachim

Hi Joachim,

Thanks for the answer.
I will implement it and see how it works.

Hi Joachim,

I was trying the HDG as you suggested but it looks like im still doing something very wrong. Now it crashes when I try to define the Facet space. Might it be a problem of 1D mesh?

The message I obtain when I try to create the Facet space is

setting all_dofs_together 0
Segmentation fault: 11
My space definition looks like
Y = H1(mesh, order=ord, dirichlet=[1,2])
F = FacetFESpace(mesh, order=ord)
and the second line is where the code collapses. Do you have an idea what might be the reason?

Best wishes and thanks again,

Hi Bart,

HDG was not fully implemented for 1D. An updated version as available on Github and you can pull it from there. If you use an installer, you can update/download ngsolve_nightly tomorrow.

If you want to use HDG, you have to use an discontinuous finite element space for your element variable (e.g. L2). For a 1d mesh, you just have one degree of freedom on the facets, thus you have a 0th order facet finite element space.

order = 4
V1 = L2(mesh, order=order)
V2 = FacetFESpace(mesh, order=0, dirichlet=[1])
V = FESpace([V1,V2])

A full example solving a laplace problem with HDG is attached.


Attachment: make_1d_mesh.py

Attachment: facetfe1d.py

Thanks for the previous answer. I defined my spaces as you advised and now they look like

Y = L2(mesh, order=ord) # space can be fixed for the whole routine Z = FacetFESpace(mesh, order=0,dirichlet = [1,2]) X = FESpace([Y,Z,Y,Z,Y,Z]) .
As always there is a but. Now I struggle with the definitions of the BC and the initial guess for the Newton loop.
What I have is

u.components[1].vec[nel]  = Vapl
u.components[3].vec[nel]=z1*beta*Vapl  + log(bc1)
u.components[5].vec[nel]=z2*beta*Vapl  + log(bc1)

which (as I believe) provides the right boundary conditions for the facet space and should provide the BC for the Newton loop. Unfortunately this way the code gets stuck on the line

a.Apply(u.vec, res)

[color=red]I fixed this issue and it is not stuck any more. I needed to increse the heapsize by seting a.Apply(u.vec, res, heapsize=100000000) Now I end up with a singular matrix so I believe the BC/initialization problem stays the same. [/color]

Is there a better way to implement the BC and the initial guess than the one I’m using?

Best wishes,


for the 1d case, your way of setting the boundary conditions should work. (The Set() function which you used before does not support points by now.)

Which matrix is singular? The result of “a.Apply(…)” is a vector.

If you do not set anything else than those three boundary values, all other values in “u.vec” are 0. Is that the vector you want to start with?



I thought that I’m doing something wrong with the Set function. Is good to know that it is not only me.

For the BC and the initial Guess.
I set 3 boundary conditions manually by


so that the facet space vector is 0 everywhere but the boundary. I don’t set anything for the values of the functions in the L2 space.

Then my system enters the Newtons loop and what I get is
linear form done
Iteration 0
<A u 0 , A u 0 >{-1}^0.5 = 190.3901437327294
Iteration 1
<A u 1 , A u 1 >
{-1}^0.5 = 2.0534197662160665

UMFPACK V5.7.4 (Feb 1, 2016): WARNING: matrix is singular

so I get a singular matrix in my Newton loop which looks like

res = u.vec.CreateVector()
du = u.vec.CreateVector()
################################# newton loop
for it in range(250):
a.Apply(u.vec, res)
du.data = a.mat.Inverse(X.FreeDofs()) * res
u.vec.data -= du
#stopping criteria
stopcritval = sqrt(abs(InnerProduct(du,res)))
print (“Iteration”,it)
print (“<A u”,it,“, A u”,it,">_{-1}^0.5 = ", stopcritval)

[/code]. The problematic line is the inverse now. Any ideas why the HDG formulation might lead to a singular matrix? I was solving the same system just in H1 space and I didn’t have this problem.

Best wishes,


could you please attach your updated file with the HDG formulation.



The problem with the singular matrices arises because in my bilinear form I have some functions that are x dependent A(x) and Sigma(x) (you can see the code attached). Before I was just assigning values to the vectors manually using Sigma.vec[i] = g(i) . As now Sigma is a GridFunction(L2) this way is a no longer valid method of assigning values. I changed it then and in the code, I attach I just set

to have something. Now the Newton loop converges but my ability to assign values is somehow limited. We are somehow back to the Set function we discussed before. Is there any way to assign the values of this functions (A and sigma) such that I can use for example indicator functions over x?

Best wishes,

Attachment: PNP_2018-01-25.py


if I understand you correct, you just want do set more complicated functions for A and Sigma?
You can use a indicator function


or any other “CoefficentFunction”.
If you call


within python, you get some information what you can do with CoefficientFunctions.


Thanks Christoph! I was looking for that for quite some time. I will try then to put it all together now.

Best wishes,

Dear Christoph,

I’ve managed to make my HDG with the Newton loop code running. As it is only a part of a bigger project I experienced some new problem. In order to use the command a.Apply where a is the bilinear form for the HDG formulation I need to increase the heap size by SetHeapSize(100000000). That works fine when I try to solve my equation once. Unfortunately, I would love to use it as a part of the som optimization problem and need to solve the equation in a loop with different parameters. When I solve the equation for the second time I encounter the problem with the heap size- the program freezes at the application of the form as it was because with the standard size of the heap. I believe I do something wrong with the memory management. I was trying to use the del and gc.collect() to free some memory but it did not help.

Any idea what might be the issue?
Best wishes,

Attachment: PNP_small.py

Hi Bart,

there was a bug in the FacetFESpace for 1d. The fix should be available in the next days.

A few comments to your code:
[li]If you want to visualize a 1d function you have to use Gridfunctions. Drawing CoefficientFunctions in 1d is not working. In case you don’t want to draw your functions (like Sigma, R, A, xpoints, …) a would advise to use the as CoefficientFunctions.Sigma = IfPos(x-1200,1,0)*IfPos(-x+2400,1,0)[/li]
defines a CoefficientFunction which you can use as parameter in your BilinearForm.
[li]When you calculate your output functions, you use J1.vec[int(nel/2)] I guess you want the function value in the middle? If so, just leave c1, gradu1 and J1 as CoefficientFunctions and evaluate at the correct point. For the domain [0,length] this can be done with J1(mesh(length/2))
The entry of the vector is just a coefficient value and might be different than the function value.
[li]You do not have to delete all the variables (you del … statement). This should be done when the PNP_solver is finished.[/li]


Dear Christoph,

I downloaded the new version of the library and now there is no problem with the heap size. So the update seems to work.

For the Grid/Coefficient function - I struggled a lot with plotting them and passing into my functions. As form time to time I want to plot some of them it is better to keep them as they are for now.

On more question - I would like to introduce some box constrain for a function and limit it value to (-2,2). I was trying with the Ifpos function but I’m missing the absolute value. What I would like to do is

Sigma. Set( IfPos(gdgd-4,2abs(gd),gd))

but the absolute value is not working. Any idea how to go over this problem?

Best wishes,

You can use the Norm CoefficientFunction, abs does not work with CFs (but I see no reason why, so I think we could add it as well…)
But if you want to limit it’s values to [-2,2] you would have to cut them off at the top and the bottom:

Sigma. Set( IfPos(gd-2,2, IfPos(2-gd,gd,-2)))


You could still use Norm as well:

Sigma.Set( IfPos(gd*gd-4,2*gd/Norm(gd),gd))

does what you want.



I was just wondering if there are any plans to implement linearization for DG? Or is support still predominantly for HDG?