--- Begin Message ---
------------------------------------------------------------
revno: 6357
committer: Johan Hake <hake.dev@xxxxxxxxx>
branch nick: work
timestamp: Mon 2011-10-17 19:04:45 -0700
message:
Fix subdomain assemble for SystemAssembler
-- Probably some possibilities for code reuse between the three Assembler files
-- It make sense to put subdomain and tabulate tensor logic into UFC
modified:
ChangeLog
dolfin/fem/SystemAssembler.cpp
test/unit/fem/python/SystemAssembler.py
--
lp:dolfin
https://code.launchpad.net/~dolfin-core/dolfin/main
Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/main/+edit-subscription
=== modified file 'ChangeLog'
--- ChangeLog 2011-10-17 23:20:50 +0000
+++ ChangeLog 2011-10-18 02:04:45 +0000
@@ -1,3 +1,4 @@
+ - Fix subdomain assemble for SystemAssembler
- Fix OpenMp assemble of scalars
- Make OpenMP assemble over sub domains work
- DirichletBC.get_boundary_values, FunctionSpace.collapse now return a dict in Python
=== modified file 'dolfin/fem/SystemAssembler.cpp'
--- dolfin/fem/SystemAssembler.cpp 2011-10-04 20:35:02 +0000
+++ dolfin/fem/SystemAssembler.cpp 2011-10-18 02:04:45 +0000
@@ -96,9 +96,33 @@
// FIXME: Some things can be simplified since we know it's a matrix and a vector
- // Sub domains not supported
- if (cell_domains || exterior_facet_domains || interior_facet_domains)
- error("SystemAssembler does not yet support subdomains.");
+ // Get cell domains
+ if (!cell_domains || cell_domains->size() == 0)
+ {
+ cell_domains = a.cell_domains_shared_ptr().get();
+ if (!cell_domains)
+ cell_domains = a.mesh().domains().cell_domains(a.mesh()).get();
+ }
+
+ // Get exterior facet domains
+ if (!exterior_facet_domains || exterior_facet_domains->size() == 0)
+ {
+ exterior_facet_domains = a.exterior_facet_domains_shared_ptr().get();
+ if (!exterior_facet_domains)
+ exterior_facet_domains = a.mesh().domains().facet_domains(a.mesh()).get();
+ }
+
+ // FIXME: Sub domains for interior facets not supported
+ // Get interior facet domains
+ //if (!interior_facet_domains || interior_facet_domains->size() == 0)
+ //{
+ // interior_facet_domains = a.interior_facet_domains_shared_ptr().get();
+ // if (!interior_facet_domains)
+ // interior_facet_domains = a.mesh().domains().facet_domains(a.mesh()).get();
+ //}
+
+ if (interior_facet_domains)
+ error("SystemAssembler does not yet support subdomains for interior facets.");
// Check arguments
AssemblerTools::check(a);
@@ -220,6 +244,7 @@
const Mesh& mesh = a.mesh();
assert(mesh.ordered());
+ const uint D = mesh.topology().dim();
// Form ranks
const uint a_rank = a.rank();
@@ -239,22 +264,53 @@
std::vector<const std::vector<uint>* > L_dofs(L_rank);
// Cell integrals
- const ufc::cell_integral* A_integral(0);
- const ufc::cell_integral* b_integral(0);
+ ufc::cell_integral* A_integral(0);
+ ufc::cell_integral* b_integral(0);
if (A_ufc.cell_integrals.size() > 0)
A_integral = A_ufc.cell_integrals[0].get();
if (b_ufc.cell_integrals.size() > 0)
b_integral = b_ufc.cell_integrals[0].get();
- // FIXME: Subdomains not supported here!
+ // Facet_integrals
+ ufc::exterior_facet_integral* A_facet_integral(0);
+ ufc::exterior_facet_integral* b_facet_integral(0);
+ if (A_ufc.exterior_facet_integrals.size() > 0)
+ A_facet_integral = A_ufc.exterior_facet_integrals[0].get();
+ if (b_ufc.exterior_facet_integrals.size() > 0)
+ b_facet_integral = b_ufc.exterior_facet_integrals[0].get();
+ // Init facet connection if we have facet integrals
+ if (A_ufc.exterior_facet_integrals.size() > 0 ||
+ b_ufc.exterior_facet_integrals.size() > 0)
+ {
+ mesh.init(D - 1);
+ mesh.init(D - 1, D);
+ }
+
Progress p("Assembling system (cell-wise)", mesh.num_cells());
for (CellIterator cell(mesh); !cell.end(); ++cell)
{
+
+ // Get integrals for sub domain (if any)
+ if (cell_domains && cell_domains->size() > 0)
+ {
+ const uint cell_domain = (*cell_domains)[*cell];
+ if (cell_domain < A_ufc.form.num_cell_domains())
+ A_integral = A_ufc.cell_integrals[cell_domain].get();
+ else
+ A_integral = 0;
+
+ if (cell_domain < b_ufc.form.num_cell_domains())
+ b_integral = b_ufc.cell_integrals[cell_domain].get();
+ else
+ b_integral = 0;
+ }
+
// Add cell tensor for A
if (A_integral)
{
- // Update to current cell
+ //info("Cell domain A: %d", (*cell_domains)[*cell]);
+ // Update to current cell
A_ufc.update(*cell);
// Tabulate cell tensor
@@ -262,10 +318,18 @@
for (uint i = 0; i < data.Ae.size(); ++i)
data.Ae[i] = A_ufc.A[i];
}
+ else
+ {
+ // Zero out tensor when no cell integral
+ for (uint i = 0; i < data.Ae.size(); ++i)
+ data.Ae[i] = 0.0;
+ }
+
// Add cell tensor for b
if (b_integral)
{
+ //info("Cell domain B: %d", (*cell_domains)[*cell]);
// Update to current cell
b_ufc.update(*cell);
@@ -274,44 +338,71 @@
for (uint i = 0; i < data.be.size(); ++i)
data.be[i] = b_ufc.A[i];
}
+ else
+ {
+ // Zero out tensor when no cell integral
+ for (uint i = 0; i < data.be.size(); ++i)
+ data.be[i] = 0.0;
+ }
// FIXME: Is assembly over facets more efficient?
// Compute exterior facet integral if present
if (A_ufc.form.num_exterior_facet_domains() > 0 ||
b_ufc.form.num_exterior_facet_domains() > 0)
{
+
+ // Get connectivity
+ const MeshConnectivity& connectivity = mesh.topology()(D, D - 1);
+
for (FacetIterator facet(*cell); !facet.end(); ++facet)
{
- // Assemble if we have an exterior facet
- if (facet->exterior())
- {
- const uint local_facet = cell->index(*facet);
- if (A_ufc.form.num_exterior_facet_domains() > 0)
- {
- A_ufc.update(*cell, local_facet);
-
- const ufc::exterior_facet_integral*
- A_facet_integral = A_ufc.exterior_facet_integrals[0].get();
- A_facet_integral->tabulate_tensor(&A_ufc.A[0],
- A_ufc.w(),
- A_ufc.cell,
- local_facet);
- for (uint i = 0; i < data.Ae.size(); i++)
- data.Ae[i] += A_ufc.A[i];
- }
- if (b_ufc.form.num_exterior_facet_domains() > 0)
- {
- b_ufc.update(*cell, local_facet);
-
- const ufc::exterior_facet_integral*
- b_facet_integral = b_ufc.exterior_facet_integrals[0].get();
- b_facet_integral->tabulate_tensor(&b_ufc.A[0],
- b_ufc.w(),
- b_ufc.cell,
- local_facet);
- for (uint i = 0; i < data.be.size(); i++)
- data.be[i] += b_ufc.A[i];
- }
+ // Only consider exterior facets
+ if (!facet->exterior())
+ continue;
+
+ const uint local_facet = cell->index(*facet);
+
+ // Get integrals for sub domain (if any)
+ if (exterior_facet_domains && exterior_facet_domains->size() > 0)
+ {
+
+ // Get global facet index
+ const uint facet_index = connectivity(cell->index())[local_facet];
+ const uint facet_domain = (*exterior_facet_domains)[facet_index];
+
+ if (facet_domain < A_ufc.form.num_exterior_facet_domains())
+ A_facet_integral = A_ufc.exterior_facet_integrals[facet_domain].get();
+ else
+ A_facet_integral = 0;
+
+ if (facet_domain < b_ufc.form.num_exterior_facet_domains())
+ b_facet_integral = b_ufc.exterior_facet_integrals[facet_domain].get();
+ else
+ b_facet_integral = 0;
+ }
+
+ // Tabualte tensor if facet integrals
+ if (A_facet_integral)
+ {
+ A_ufc.update(*cell, local_facet);
+ A_facet_integral->tabulate_tensor(&A_ufc.A[0],
+ A_ufc.w(),
+ A_ufc.cell,
+ local_facet);
+ for (uint i = 0; i < data.Ae.size(); i++)
+ data.Ae[i] += A_ufc.A[i];
+ }
+
+ if (b_facet_integral)
+ {
+ b_ufc.update(*cell, local_facet);
+ b_facet_integral->tabulate_tensor(&b_ufc.A[0],
+ b_ufc.w(),
+ b_ufc.cell,
+ local_facet);
+ for (uint i = 0; i < data.be.size(); i++)
+ data.be[i] += b_ufc.A[i];
+
}
}
}
=== modified file 'test/unit/fem/python/SystemAssembler.py'
--- test/unit/fem/python/SystemAssembler.py 2011-10-05 11:14:14 +0000
+++ test/unit/fem/python/SystemAssembler.py 2011-10-18 02:04:45 +0000
@@ -96,7 +96,7 @@
"Test assembly over subdomains with markers stored as part of mesh"
# Create a mesh of the unit cube
- mesh = UnitCube(4, 4, 4)
+ mesh = UnitCube(2, 2, 2)
# Define subdomains for 3 faces of the unit cube
class F0(SubDomain):
@@ -137,23 +137,28 @@
a1 = 1*u*v*dx(1) + 2*u*v*ds(3)
L1 = 1*v*dx(1) + 2*v*ds(3)
- # FIXME: Currently disabled, need to fix bug
-
# Used for computing reference values
- #A0 = assemble(a0)
- #b0 = assemble(L0)
- #A1 = assemble(a1)
- #b1 = assemble(L1)
+ A0 = assemble(a0)
+ b0 = assemble(L0)
+ A1 = assemble(a1)
+ b1 = assemble(L1)
+
+ A0_ref = A0.norm("frobenius")
+ A1_ref = A1.norm("frobenius")
+ b0_ref = b0.norm("l2")
+ b1_ref = b1.norm("l2")
# Assemble system
- #A0, b0 = assemble_system(a0, L0)
- #A1, b1 = assemble_system(a1, L1)
+ A0, b0 = assemble_system(a0, L0)
+ A1, b1 = assemble_system(a1, L1)
- # Assemble and check values
- #self.assertAlmostEqual(A0.norm("frobenius"), 0.693043954566, 10)
- #self.assertAlmostEqual(b0.norm("l2"), 1.28061997552, 10)
- #self.assertAlmostEqual(A1.norm("frobenius"), 0.45406526606, 10)
- #self.assertAlmostEqual(b1.norm("l2"), 0.84277689513, 10)
+ # FIXME: Should these test work? (They don't...) The matrix and the
+ # FIXME: vector from AssembleSystem need to be different from that
+ # FIXME: of Assemble
+ self.assertAlmostEqual(A0.norm("frobenius"), A0_ref, 10)
+ self.assertAlmostEqual(A1.norm("frobenius"), A1_ref, 10)
+ self.assertAlmostEqual(b0.norm("l2"), b0_ref, 10)
+ self.assertAlmostEqual(b1.norm("l2"), b1_ref, 10)
if __name__ == "__main__":
print ""
--- End Message ---
