Commits

Anders Logg committed 92b4889

More work on CCFEM sparsity pattern building and assembly. Basic
functionality in place (but still very incomplete).

Comments (0)

Files changed (11)

demo/undocumented/ccfem-poisson/cpp/main.cpp

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-06-26
-// Last changed: 2013-09-18
+// Last changed: 2013-09-24
 //
 // This demo program solves Poisson's equation using a Cut and
 // Composite Finite Element Method (CCFEM) on a domain defined by
   info("THIS DEMO IS WORK IN PROGRESS!");
 
   // Increase log level
-  set_log_level(PROGRESS);
+  set_log_level(DBG);
 
   // Create meshes
-  UnitSquareMesh square(8, 8);
-  UnitCircleMesh circle_1(8);
-  UnitCircleMesh circle_2(8);
+  UnitSquareMesh square(4, 4);
+  UnitCircleMesh circle_1(3);
+  UnitCircleMesh circle_2(3);
 
   // Displace circle meshes
   Point dx(0.5, 0.5);
   CCFEMAssembler assembler;
   assembler.assemble(A, a);
 
+  info(A, true);
+
   return 0;
 }

dolfin/fem/Assembler.cpp

 // Modified by Martin Alnaes 2013
 //
 // First added:  2007-01-17
-// Last changed: 2013-01-18
+// Last changed: 2013-09-19
 
 #include <boost/scoped_ptr.hpp>
 

dolfin/fem/Assembler.h

 // Modified by Joachim B Haga 2012
 //
 // First added:  2007-01-17
-// Last changed: 2013-09-12
+// Last changed: 2013-09-19
 
 #ifndef __ASSEMBLER_H
 #define __ASSEMBLER_H

dolfin/fem/AssemblerBase.cpp

 // Modified by Johannes Ring, 2012
 //
 // First added:  2007-01-17
-// Last changed: 2013-09-18
+// Last changed: 2013-09-19
 
 #include <boost/scoped_ptr.hpp>
 #include <dolfin/common/Timer.h>
 #include "SparsityPatternBuilder.h"
 #include "AssemblerBase.h"
 
-
 using namespace dolfin;
 
 //-----------------------------------------------------------------------------

dolfin/fem/AssemblerBase.h

 // Modified by Ola Skavhaug, 2008.
 //
 // First added:  2007-01-17
-// Last changed: 2009-10-06
+// Last changed: 2013-09-19
 
 #ifndef __ASSEMBLER_BASE_H
 #define __ASSEMBLER_BASE_H
   public:
 
     // Check form
-    AssemblerBase() : reset_sparsity(true), add_values(false),
-        finalize_tensor(true), keep_diagonal(false) {}
+    AssemblerBase() :
+      reset_sparsity(true),
+      add_values(false),
+      finalize_tensor(true),
+      keep_diagonal(false) {}
 
     /// reset_sparsity (bool)
     ///     Default value is true.

dolfin/fem/CCFEMAssembler.cpp

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-09-12
-// Last changed: 2013-09-19
+// Last changed: 2013-09-24
 
 #include <dolfin/log/log.h>
 #include <dolfin/la/GenericTensor.h>
+#include <dolfin/la/GenericMatrix.h>
 #include <dolfin/la/GenericLinearAlgebraFactory.h>
 #include <dolfin/la/TensorLayout.h>
 #include <dolfin/function/CCFEMFunctionSpace.h>
 
 #include "SparsityPatternBuilder.h"
+#include "UFC.h"
+#include "Form.h"
 #include "CCFEMForm.h"
+#include "CCFEMDofMap.h"
 #include "CCFEMAssembler.h"
 
 using namespace dolfin;
 
 //-----------------------------------------------------------------------------
+CCFEMAssembler::CCFEMAssembler()
+{
+  // Do nothing
+}
+//-----------------------------------------------------------------------------
 void CCFEMAssembler::assemble(GenericTensor& A, const CCFEMForm& a)
 {
   begin(PROGRESS, "Assembling tensor over CCFEM function space.");
   // Initialize global tensor
   init_global_tensor(A, a);
 
+  // Assemble over cells
+  assemble_cells(A, a);
+
+  // Finalize assembly of global tensor
+  A.apply("add");
+
   end();
 }
 //-----------------------------------------------------------------------------
+void CCFEMAssembler::assemble_cells(GenericTensor& A, const CCFEMForm& a)
+{
+  // Note: This implementation does not yet handle subdomains.
+
+  log(PROGRESS, "Assembling CCFEM form over cells.");
+
+  // Form rank
+  const std::size_t form_rank = a.rank();
+
+  // Vector to hold dof map for a cell
+  std::vector<const std::vector<dolfin::la_index>* > dofs(form_rank);
+
+  // Collect pointers to dof maps
+  std::vector<const CCFEMDofMap*> dofmaps;
+  for (std::size_t i = 0; i < form_rank; i++)
+    dofmaps.push_back(a.function_space(i)->dofmap().get());
+
+  // Iterate over parts
+  for (std::size_t part = 0; part < a.num_parts(); part++)
+  {
+    // Get form for current part
+    const Form& a_part = *a.part(part);
+
+    // Set current part for dofmaps
+    for (std::size_t i = 0; i < form_rank; i++)
+      dofmaps[i]->set_current_part(part);
+
+    // Create data structure for local assembly data
+    UFC ufc(a_part);
+
+    // Extract mesh
+    const Mesh& mesh = a_part.mesh();
+
+    // Skip assembly if there are no cell integrals
+    if (!ufc.form.has_cell_integrals())
+      return;
+
+    // Cell integral
+    ufc::cell_integral* integral = ufc.default_cell_integral.get();
+
+    // Iterate over cells
+    for (CellIterator cell(mesh); !cell.end(); ++cell)
+    {
+      // Update to current cell
+      ufc.update(*cell);
+
+      // Get local-to-global dof maps for cell
+      for (std::size_t i = 0; i < form_rank; ++i)
+        dofs[i] = &(dofmaps[i]->cell_dofs(cell->index()));
+
+      // Tabulate cell tensor
+      integral->tabulate_tensor(&ufc.A[0], ufc.w(),
+                                &ufc.cell.vertex_coordinates[0],
+                                ufc.cell.orientation);
+
+      // Add entries to global tensor
+      A.add(&ufc.A[0], dofs);
+    }
+  }
+}
+//-----------------------------------------------------------------------------
 void CCFEMAssembler::init_global_tensor(GenericTensor& A, const CCFEMForm& a)
 {
   log(PROGRESS, "Initializing global tensor.");
   {
     GenericSparsityPattern& pattern = *tensor_layout->sparsity_pattern();
     SparsityPatternBuilder::build_ccfem(pattern, a);
+
+    // FIXME: debugging
+    cout << tensor_layout->str() << endl;
   }
+
+  // Initialize tensor
+  A.init(*tensor_layout);
+
+  // Insert zeros on the diagonal as diagonal entries may be prematurely
+  // optimised away by the linear algebra backend when calling
+  // GenericMatrix::apply, e.g. PETSc does this then errors when matrices
+  // have no diagonal entry inserted.
+  if (A.rank() == 2)
+  {
+    // Down cast to GenericMatrix
+    GenericMatrix& _A = A.down_cast<GenericMatrix>();
+
+    // Loop over rows and insert 0.0 on the diagonal
+    const double block = 0.0;
+    const std::pair<std::size_t, std::size_t> row_range = A.local_range(0);
+    const std::size_t range = std::min(row_range.second, A.size(1));
+    for (std::size_t i = row_range.first; i < range; i++)
+    {
+      dolfin::la_index _i = i;
+      _A.set(&block, 1, &_i, 1, &_i);
+    }
+    A.apply("flush");
+  }
+
+  // Set tensor to zero
+  A.zero();
 }
 //-----------------------------------------------------------------------------

dolfin/fem/CCFEMAssembler.h

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // First added:  2013-09-12
-// Last changed: 2013-09-12
+// Last changed: 2013-09-24
 
 #ifndef __CCFEM_ASSEMBLER_H
 #define __CCFEM_ASSEMBLER_H
 
 #include "AssemblerBase.h"
+#include "Assembler.h"
 
 namespace dolfin
 {
   public:
 
     /// Constructor
-    CCFEMAssembler() {}
+    CCFEMAssembler();
 
     /// Assemble tensor from given form
     ///
 
   private:
 
+    // Assemble over cells
+    void assemble_cells(GenericTensor& A, const CCFEMForm& a);
+
     // Initialize global tensor
     void init_global_tensor(GenericTensor& A, const CCFEMForm& a);
 

dolfin/fem/SparsityPatternBuilder.cpp

 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // Modified by Ola Skavhaug 2007
-// Modified by Anders Logg 2008-2011
+// Modified by Anders Logg 2008-2013
 //
 // First added:  2007-05-24
-// Last changed: 2013-09-19
+// Last changed: 2013-09-24
 
 #include <dolfin/common/timing.h>
 #include <dolfin/common/MPI.h>
                                    bool cells,
                                    bool interior_facets,
                                    bool exterior_facets,
-                                   bool diagonal)
+                                   bool diagonal,
+                                   bool init,
+                                   bool finalize)
 {
   const std::size_t rank = dofmaps.size();
 
   }
 
   // Initialise sparsity pattern
-  sparsity_pattern.init(global_dimensions, local_range, off_process_owner);
+  if (init)
+    sparsity_pattern.init(global_dimensions, local_range, off_process_owner);
 
   // Only build for rank >= 2 (matrices and higher order tensors) that
   // require sparsity details
   }
 
   // Finalize sparsity pattern (communicate off-process terms)
-  sparsity_pattern.apply();
+  if (finalize)
+    sparsity_pattern.apply();
 }
 //-----------------------------------------------------------------------------
 void SparsityPatternBuilder::build_ccfem(GenericSparsityPattern& sparsity_pattern,
     // Get mesh on current part (assume it's the same for all arguments)
     const Mesh& mesh = *form.function_space(0)->part(part)->mesh();
 
+    // Check whether to initialize or finalize sparsity pattern
+    const bool init = part == 0;
+    const bool finalize = part == form.num_parts() - 1;
+
     // Build sparsity pattern for part
-    build(sparsity_pattern, mesh, dofmaps, true, false, false, true);
+    build(sparsity_pattern, mesh, dofmaps, true, false, false, true,
+          init, finalize);
   }
 }
 //-----------------------------------------------------------------------------

dolfin/fem/SparsityPatternBuilder.h

 // Modified by Anders Logg 2008-2013
 //
 // First added:  2007-05-24
-// Last changed: 2013-09-18
+// Last changed: 2013-09-24
 
 #ifndef __SPARSITY_PATTERN_BUILDER_H
 #define __SPARSITY_PATTERN_BUILDER_H
                       bool cells,
                       bool interior_facets,
                       bool exterior_facets,
-                      bool diagonal);
+                      bool diagonal,
+                      bool init=true,
+                      bool finalize=true);
 
     /// Build sparsity pattern for assembly of given CCFEM form
     static void build_ccfem(GenericSparsityPattern& sparsity_pattern,

dolfin/la/PETScMatrix.cpp

 // Modified by Jan Blechta 2013
 //
 // First added:  2004
-// Last changed: 2012-03-15
+// Last changed: 2013-09-24
 
 #ifdef HAS_PETSC
 
     std::vector<std::size_t> num_nonzeros(M);
     sparsity_pattern.num_nonzeros_diagonal(num_nonzeros);
 
+    // FIXME: Debugging
+    //cout << "Number of nonzeros for PETSc matrix:" << endl;
+    //for (unsigned int i = 0; i < M; i++)
+    //  cout << "i = " << i << ": " << num_nonzeros[i] << endl;
+
     // Create matrix
     ierr = MatCreate(PETSC_COMM_SELF, _A.get());
     if (ierr != 0) petsc_error(ierr, __FILE__, "MatCreate");
 // Modified by Benjamin Kehlet, 2012
 //
 // First added:  2007-05-02
-// Last changed: 2012-09-16
+// Last changed: 2013-09-19
 
 #include <cstdlib>
 #include <sstream>
     wrapper(new ExpressionWrapper(expression, mesh));
   return plot_object(wrapper, p, VTKPlotter::to_key(*expression));
 }
+//-----------------------------------------------------------------------------