Dirichlet boundary - enforce(A, b, D=D)

[JarosT55]

Hi,

I have looked at the behaviour of $enforce(A, b, D=D)$.
It takes a symmetric array $A( A(i,j) \equiv A(j,i) )$ and set

$\quad A(i,j)=A(j,i)=0.0 \quad \forall \implies dof(i,j) ∈ D \land i \neq j.$
$\quad A(i,j)=A(j,i)=1.0 \quad \forall \implies dof(i,j) ∈ D \land i = j.$

To my surprise, the procedure $enforce(A, b, D=D)# enforces only 'rows' of A.
Thus, as a consequence, $A$ loses its symmetricity.

In most case $A$ is a subject of conversion to $C = condense(A, b, D=D)$,
which once again produces $C$ as a symmetric array ($C[0]$).

Does $enforce(A, b, D=D)$ in its current form has other function/purpose ... I can not comprehend?

As per the function doc label:
"A linear system with the enforced rows/diagonals set to zero/one."

I reiterate, ... why only rows and not columns also?

Regards
MichaelT

[kinnala]

This does not and should not lead to a symmetric matrix.

Edit: Setting also columns to zero will lead to incorrect results.

[JarosT55]

Hi Mr Kinnala,

Thank you for the input.

I think I know the reason now. Doing some mackageegee around Dirichlet conditions, I was fixated on $A$ matrix, and the loss of its symmetricity after the $enforce()$ has concerned me.
Now I understand the purpose of $enforce()$. It is a way of preparing linear system $A \cdot b$ for a successive step of solving it and/or $LU$ decomposition of the linear system while $D$ dofs have been removed.

Regards
MichaelT

[JarosT55]

Hi Mr Kinnala,

I have done some performance testing of:
$A, b=enforce(A,b,D)$
$x=solve(A,b)$
against:
$x=solve(*condense(A,b,D))$

The results ($x$) are identical. The difference $solve(A,_)$ in $A$ I see is LU vs symmetric matrix.
Both, the calculation timings and memory use are statistically similar.
Although there might be some uses for $A$ in $LU$ form ... for me a symmetric matrix representing mesh's neighbourhood is more natural in a physical sense, as path between connected nodes is bidirectional with equal 'speeds' in both ways.
The situation might be different when, e.g. space/time graph (mesh) is involved, then $LU$ arrays might come to the rescue.
However, we have a 'normal' situation when 1 hour in the past equals 1 hour in the future, even while traveling fast through our nonrelativistic space.

Test results:

Duration (A.shape=(263169, 263169), loop=10):
$A, b=enforce(A,b,D)$
$x1=solve(A,b)$
Time: 0:05:04.639892 MemoryPeak 44,519,205

Duration (A.shape=(263169, 263169), loop=10):
$x2=solve(*condense(A,b,D))$
Time: 0:04:32.623942 MemoryPeak 40,787,255

Results equality:
$(x1 EQ x2) = True$

Test code, based on ex01.

# ex01

from skfem import *
from skfem.helpers import dot, grad
from datetime import datetime
import tracemalloc
import gc

# # enable additional mesh validity checks, sacrificing performance
# import logging
# logging.basicConfig(format='%(levelname)s %(asctime)s %(name)s %(message)s')
# logging.getLogger('skfem').setLevel(logging.DEBUG)

# create the mesh
m = MeshTri().refined(6)
# or, with your own points and cells:
# m = MeshTri(points, cells)

e = ElementTriP1()
basis = Basis(m, e)

# this method could also be imported from skfem.models.laplace
@BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))


# this method could also be imported from skfem.models.unit_load
@LinearForm
def rhs(v, _):
    return 1.0 * v

def visualize():
    from skfem.visuals.matplotlib import plot
    return plot(m, x, shading='gouraud', colorbar=True)

A = asm(laplace, basis)
b = asm(rhs, basis)
# or:
# A = laplace.assemble(basis)
# b = rhs.assemble(basis)

# enforce Dirichlet boundary conditions
A, b = enforce(A, b, D=m.boundary_nodes())
# solve -- can be anything that takes a sparse matrix and a right-hand side
x = solve(A, b)
visualize().show()

# The equivalent version without enforce()  
# solve -- can be anything that takes a sparse matrix and a right-hand side
x = solve(*condense(A, b, D=m.boundary_nodes()))
visualize().show()    

# Do compare code execution timing
re= 8            # refined(8) ~ 20 sec
m = MeshTri().refined(re)
e = ElementTriP1()
D = m.boundary_nodes()
basis = Basis(m, e)
A = asm(laplace, basis)
b = asm(rhs, basis)
nl = 10
gc.collect()
#--------------------------------------------------
start_time = datetime.now()
gc.collect()
tracemalloc.start()
for i in range(nl):
    A, b = enforce(A, b, D=D)
    x1 = solve(A, b)
    del(x1)
    gc.collect()
memoryTrace = tracemalloc.get_traced_memory()
tracemalloc.stop()    
end_time = datetime.now()
print('Duration (A.shape={}, loop={}):\n    A, b=enforce(A,b,D)\n    x=solve(A,b)\nTime: {}   MemoryPeak {:,}\n'.format(A.shape, nl, end_time - start_time, memoryTrace[1]))
#--------------------------------------------------
start_time = datetime.now()
gc.collect()
tracemalloc.start()
for i in range(nl):
    x2 = solve(*condense(A, b, D=D))
    del(x2)
    gc.collect()
memoryTrace = tracemalloc.get_traced_memory()
tracemalloc.stop()    
end_time = datetime.now()
print('Duration (A.shape={}, loop={}):\n    x=solve(*condense(A,b,D))\nTime: {}   MemoryPeak {:,}\n'.format(A.shape, nl, end_time - start_time, memoryTrace[1]))
#--------------------------------------------------
print('(x1 EQ x2) =', str(np.allclose(x1, x2, rtol=1e-05, atol=1e-08))) 

`

Regards
MichaelT