# Symmetric DirichletBC?

**URL:** https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659
**Category:** Uncategorized
**Created:** [April 23, 2019, 5:05pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659 "2019-04-23T17:05:25Z")
**Posts on this page:** 11
**Page:** 1

<div class="post-metadata">

### Author: ![nice-butt](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nice-butt/32/161_2.png) [@nice-butt](https://fenicsproject.discourse.group/u/nice-butt)
#### Post date: [April 23, 2019, 5:05pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/1 "2019-04-23T17:05:25Z")

</div>

Hi

I am solving nonlinear solid mechanics computations. So I overrides `dolfin.NonlinearProblem` to apply assembly functionalities as well as applying DirichletBC using `bc.apply`

When I check if the system Jacobian matrix is symmetric (`A.mat().isSymmetric(...)`), it says no. Yet the formulation tells the symmetry. It happens to the case of solving nonlinear Poisson’s equation.

Does anyone have an idea how to get DirichletBC applied without breaking symmetry? Or is it just me missing something here?

Thanks,  
Victor

---

<div class="post-metadata">

### Author: ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)
#### Post date: [April 23, 2019, 6:51pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/2 "2019-04-23T18:51:35Z")

</div>

One way is to use `assemble_system`, as in the following example:

```python
from dolfin import *
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh,"CG",1)
bc = DirichletBC(V,Constant(0.0),"on_boundary")
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u),grad(v))*dx
L = Constant(1.0)*v*dx

# The `assemble_system` function maintains symmetry:
A,B = assemble_system(a,L,bcs=[bc])
print(as_backend_type(A).mat().isSymmetric())

# Using `assemble` followed by `apply` does not:
A = assemble(a)
bc.apply(A)
print(as_backend_type(A).mat().isSymmetric())

```

---

<div class="post-metadata">

### Author: ![nice-butt](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nice-butt/32/161_2.png) [@nice-butt](https://fenicsproject.discourse.group/u/nice-butt)
#### Post date: [April 23, 2019, 8:49pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/3 "2019-04-23T20:49:25Z")

</div>

Thanks, will try this out to see how i can incorporate it in subclass of NonlinearProblem

---

<div class="post-metadata">

### Author: ![nice-butt](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nice-butt/32/161_2.png) [@nice-butt](https://fenicsproject.discourse.group/u/nice-butt)
#### Post date: [April 23, 2019, 10:54pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/4 "2019-04-23T22:54:25Z")

</div>

Didn’t get any luck when using `assemble_system` along with `NonlinearProblem` for nonlinear problems. A very sloppy trial for idea testing is:

```python
class MyNP(NonlinearProblem):
    def __init__ (self, f, j, bcs):
        self.A = PETScMatrix()
        self.b = PETScVector()
        self.f = f
        self.j = j
        self.bcs = bcs
    
    def F(self, b, x):
        assemble_system(A_form=self.j, b_form=self.f, bcs=self.bcs, A_tensor=self.A, b_tensor=b)
    
    def J(self, A, x):
        assemble_system(A_form=self.j, b_form=self.f, bcs=self.bcs, A_tensor=A, b_tensor=self.b)

```

Unfortunately, it doesn’t work that residual just goes up iteration by iteration until `nan` come up and stops the simulation.

Any thing wrong?  
Thanks  
Victor

---

<div class="post-metadata">

### Author: ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)
#### Post date: [April 24, 2019, 12:19pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/5 "2019-04-24T12:19:09Z")

</div>

You can use the optional `x0` argument of `assemble_system` to set BCs in a way consistent with the system for a Newton step:

```python
from dolfin import *

mesh = UnitSquareMesh(100,100)
V = FunctionSpace(mesh,"Lagrange",1)
u = Function(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
f = x[0]*x[1]
res = (1.0+u*u)*inner(grad(u),grad(v))*dx - inner(f,v)*dx
Dres = derivative(res,u)

expr = Expression("sin(10.0*x[0]) + x[1]*x[0]",degree=3)
bc = DirichletBC(V,expr,"on_boundary")

class CustomNonlinearProblem(NonlinearProblem):
    def form(self,A,P,b,x):
        assemble_system(Dres,res,A_tensor=A,b_tensor=b,bcs=[bc,],x0=x)
    def F(self,b,x):
        pass
    def J(self,A,x):
        pass

problem = CustomNonlinearProblem()
solver = PETScSNESSolver()
solver.solve(problem,u.vector())

from matplotlib import pyplot as plt
plot(u)
plt.show()

```

---

<div class="post-metadata">

### Author: ![nice-butt](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nice-butt/32/161_2.png) [@nice-butt](https://fenicsproject.discourse.group/u/nice-butt)
#### Post date: [April 24, 2019, 4:57pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/7 "2019-04-24T16:57:56Z")

</div>

After clearing the misunderstanding of the symbol `x` that you used at different spots for different meanings, this works as a charm!

Thanks!  
Victor

---

<div class="post-metadata">

### Author: ![nice-butt](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nice-butt/32/161_2.png) [@nice-butt](https://fenicsproject.discourse.group/u/nice-butt)
#### Post date: [April 28, 2019, 11:43pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/8 "2019-04-28T23:43:04Z")

</div>

One question tho: when using the `CustomNonlinearProblem` along with `dolfin.PETScSNESSolver`, it seems to be called multiple times. To see that, I modified `form`:

```auto
def form(self, A, P, b, x):
    print('assembly starts')
    assemble_system(...)
    print('assembly ends')

```

I saw multiple times of these two lines coming out when using `PETScSNESSolver`, but not with `NewtonSolver`. I suspect that assembly of jacobian then happens multiple times per iteration, you know a resolution to avoid this?

Thanks  
Victor

---

<div class="post-metadata">

### Author: ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)
#### Post date: [April 29, 2019, 1:35pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/9 "2019-04-29T13:35:22Z")

</div>

Hmm… Looking into the code here

[https://bitbucket.org/fenics-project/dolfin/src/master/dolfin/nls/](https://bitbucket.org/fenics-project/dolfin/src/master/dolfin/nls/)

it’s clear that `form` is being called before each of `J` and `F` (and not just once, before both), although I don’t see any obvious workaround. (Trying to control calling of `assemble_system` using a Boolean that is switched off by `J` and on again by `F` doesn’t seem to work.)

---

<div class="post-metadata">

### Author: ![nice-butt](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nice-butt/32/161_2.png) [@nice-butt](https://fenicsproject.discourse.group/u/nice-butt)
#### Post date: [April 29, 2019, 4:56pm UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/10 "2019-04-29T16:56:01Z")

</div>

What I was thinking was for the first time `form` being called, I can assemble both `A` and `b`.

Case 1  
If the first call is for A, then I assemble with self.b, and when calling for assembling b, I can directly assign b with self.b.

Case 2  
If the first time is for b, then I assemble with self.A and the second call, I assign A with self.A. The reason is assembly of Jacobian for my case seems not to be trivial.

But I don’t see a way of doing for the Case 2. If this idea sounds reasonable, do you know how?

Thanks  
Victor

---

<div class="post-metadata">

### Author: ![Daniel\_Korchinski](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/daniel_korchinski/32/780_2.png) [@Daniel\_Korchinski](https://fenicsproject.discourse.group/u/Daniel_Korchinski)
#### Post date: [February 4, 2020, 2:38am UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/11 "2020-02-04T02:38:58Z")

</div>

Hi Victor,

I’m not sure that you’re still working on this problem, but something that may be of use to you (or future readers), when trying to avoid reassembly of the A matrix, is to switch from using the `assemble_system` function, and instead use a `SystemAssembler`. The SystemAssembler class is designed to apply boundary conditions at assembly (thus maintaining the symmetry in A), but can assemble b or A independently or simultaneously. You initialize an assembler as:

```auto
assembler = SystemAssembler(Dres,res,bcs=[bc,]) 
```

Then, if you want to assemble both A and b simultaneously, you can use:

```auto
assembler.assemble(A,b,x0)
```

or if you just want to reassemble the right-hand side (without having to re-assemble A):

```auto
assembler.assemble(b,x0)
```

Cheers,  
Daniel

---

<div class="post-metadata">

### Author: ![nice-butt](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nice-butt/32/161_2.png) [@nice-butt](https://fenicsproject.discourse.group/u/nice-butt)
#### Post date: [February 4, 2020, 3:01am UTC](https://fenicsproject.discourse.group/t/symmetric-dirichletbc/659/12 "2020-02-04T03:01:49Z")

</div>

Hi Daniel

Ah, I c. Will try when having a chance. Thanks!
