simple example on how to use dolfin::assemble_local


842
views
3
13 months ago by
AJ 
I would like to call dolfin FEM assemble routine for each element by providing 2d-elements and triangles (extracted from GMSH).

Does anyone have tried it and if yes can you provide a simple example code which can help me get started with it?

Thanks,
- AJ
Do you want to assemble the local matrix for each triangle, or the global matrix?  Minimal code?
written 13 months ago by pf4d 
Hello Evan,

I want to see if it is possible to call FEniCS assembly routines for each element and then build the global matrix from those local assembly calls. 

Thanks
written 13 months ago by AJ 
Interesting...  let me see....
written 13 months ago by pf4d 

2 Answers


6
13 months ago by

Here is a little code snippet which assembles a global tensors from local cell tensors using the assemble_local function.
It follows the logic of Alg. 2 in Ch. 6 of "Automated Solution of Differential Equations by the Finite Element Method".
The result is verified against the default fenics implementation by direct comparision.

# Define model problem
import numpy as np
import dolfin as df

mesh = df.RectangleMesh(df.Point(0,0), df.Point(1,1), 4, 4)
V = df.FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 2)
u = df.TrialFunction(V)
v = df.TestFunction(V)

form = df.inner(df.grad(u), df.grad(v)) * df.dx

# Assemble global matrix from local matrices
A = np.zeros([V.dim(), V.dim()])

for cell in df.cells(V.mesh()):
    A_local = df.assemble_local(form, cell)
    local_to_global_map = V.dofmap().cell_dofs(cell.index())
    
    # Copy local matrix into global matrix
    for (idx_loc_1, idx_glob_1) in enumerate(local_to_global_map):
        for (idx_loc_2, idx_glob_2) in enumerate(local_to_global_map):
            A[idx_glob_1, idx_glob_2] += A_local[idx_loc_1, idx_loc_2]

# Verification against fenics:
A_df = df.assemble(form)
A_dense = df.as_backend_type(A_df).mat().convert('dense')
assert np.allclose(A_dense.getDenseArray() - A, 0.), "Both matrices should be the same"


Hope this helps and leads you on the right track!

Oh, very nice indeed!
written 13 months ago by pf4d 
Thanks, Wagnandr. That's is working great!

Following that, how can we apply boundary conditions to the Global matrix 'A' after assembling the way you explained? 

I tried "bc.apply(A)" but get the following error:
TypeError: in method 'DirichletBC_apply', argument 2 of type 'dolfin::GenericVector &'

It works if directly applied on "A_direct = df.assemble(form)". 

Appreciate your help.
-A
       


written 11 months ago by AJ 
Hi wagnandr,

Any suggestion on how to apply Boundary Conditions on global assembly matrix (A) assembled using above method?

Thanks,
AJ
written 9 months ago by AJ 
0
13 months ago by
pf4d 
With help from old docs  :

https://fenicsproject.org/qa/927/identify-and-assemble-on-boundary-elements

from fenics import *

mesh = UnitSquareMesh(1,1)
P1   = FiniteElement('CG', triangle, 1)
Q    = FunctionSpace(mesh, P1)

# create a cell function for each cell :
cf   = CellFunction('size_t',  mesh, 0)

# for example, the first cell :
one_cell     = cells(mesh).next();
cf[one_cell] = 1

# re-define the measure :
dx = Measure('dx', subdomain_data=cf)

# variational form :
u = TestFunction(Q)
v = TrialFunction(Q)

a = inner(grad(u), grad(v))*dx(1)

# assemble over one cell :
A = assemble(a)

print A.array()
​


This gives a global matrix with non-zero values for the component of the local matrix only.

Thanks! That works. But in your example, the global mesh have just one cell. 
What if I have mesh with N cells and I want to assemble local matrix on each cell to form global matrix?

Any suggestion?
written 13 months ago by AJ 
There are two cells in this example, one upper and one lower triangle partition of the square.

I am uncertain of the best way to extract the individual local tensors; however, you can construct the global tensor by addition over all single-element assembled tensors, if that helps :

from fenics import *
import numpy as np

mesh = UnitSquareMesh(1,1)
P1   = FiniteElement('CG', triangle, 1)
Q    = FunctionSpace(mesh, P1)

# create a cell function for each cell :
cf   = CellFunction('size_t',  mesh, 0)

# iterate over all cells :
for i,cell in enumerate(cells(mesh)):
  cf[cell] = i;

# re-define the measure :
dx = Measure('dx', subdomain_data=cf)

# variational form :
u = TestFunction(Q)
v = TrialFunction(Q)

# local tensor list :
A_a = []

# assemble global contribution of local tensors :
print "LOCAL STIFFNESS:\n"
for i in range(mesh.num_cells()):
  A_i = assemble(inner(grad(u), grad(v))*dx(i))
  print A_i.array(), "\n"
  A_m = as_backend_type(A_i).mat()     # add local petsc matrix
  A_a.append(A_m)

A = np.sum(np.array(A_a))

# convert back to dolfin matrix :
C = Matrix(PETScMatrix(A))

print "GLOBAL STIFFNESS:\n", C.array()
​


Thanks : https://fenicsproject.org/qa/3150/how-to-perform-a-matrix-multiply
written 13 months ago by pf4d 
Please login to add an answer/comment or follow this question.

Similar posts:
Search »