simple example on how to use dolfin::assemble_local


589
views
3
9 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
Community: FEniCS Project
Do you want to assemble the local matrix for each triangle, or the global matrix?  Minimal code?
written 9 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 9 months ago by AJ  
Interesting...  let me see....
written 9 months ago by pf4d  

2 Answers


6
9 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 9 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 7 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 5 months ago by AJ  
0
9 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 9 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 9 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »