]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
the Simple classes now use constraints as well; some small modifications
authorjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 28 Mar 2010 17:10:52 +0000 (17:10 +0000)
committerjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 28 Mar 2010 17:10:52 +0000 (17:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@20898 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/mesh_worker_assembler.h
deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h
deal.II/deal.II/include/numerics/newton.h
deal.II/deal.II/include/numerics/newton.templates.h
deal.II/deal.II/include/numerics/operator.h
deal.II/deal.II/include/numerics/operator.templates.h
deal.II/deal.II/source/numerics/data_out.cc
deal.II/deal.II/source/numerics/operator.all_dimensions.cc
deal.II/deal.II/source/numerics/operator.all_dimensions.inst.in [new file with mode: 0644]
deal.II/deal.II/source/numerics/operator.cc [new file with mode: 0644]
deal.II/deal.II/source/numerics/operator.inst.in

index abd9c729e118789baeea3dc0c6333f2573579ced..4ec6764d872ab734447e13524d6469a270baf6c8 100644 (file)
@@ -265,6 +265,10 @@ namespace MeshWorker
     {
       public:
        void initialize(NamedData<VECTOR*>& results);
+                                        /**
+                                         * Initialize the constraints. 
+                                         */
+        void initialize(const ConstraintMatrix& constraints);
                                         /**
                                          * Initialize the local data
                                          * in the
@@ -306,6 +310,10 @@ namespace MeshWorker
                                          * filled by assemble().
                                          */
        NamedData<SmartPointer<VECTOR,ResidualSimple<VECTOR> > > residuals;
+      /**
+       * A pointer to the object containing constraints.
+       */
+      SmartPointer<const ConstraintMatrix,ResidualSimple<VECTOR> > constraints;
     };
     
 /**
@@ -394,16 +402,19 @@ namespace MeshWorker
                                          * stored as a vector of
                                          * pointers.
                                          */
-       NamedData<SmartPointer<VECTOR,ResidualLocalBlocksToGlobalBlocks<VECTOR> > > residuals;
+       NamedData<SmartPointer<VECTOR,
+          ResidualLocalBlocksToGlobalBlocks<VECTOR> > > residuals;
        
       /**
        * A pointer to the object containing the block structure.
        */
-      SmartPointer<const BlockInfo> block_info;
+      SmartPointer<const BlockInfo,
+        ResidualLocalBlocksToGlobalBlocks<VECTOR> > block_info;
       /**
        * A pointer to the object containing constraints.
        */
-      SmartPointer<const ConstraintMatrix, VECTOR> constraints;
+      SmartPointer<const ConstraintMatrix,
+        ResidualLocalBlocksToGlobalBlocks<VECTOR> > constraints;
    };
 
 
@@ -440,6 +451,10 @@ namespace MeshWorker
                                          * for later assembling.
                                          */
        void initialize(MATRIX& m);
+                                        /**
+                                         * Initialize the constraints. 
+                                         */
+        void initialize(const ConstraintMatrix& constraints);
                                         /**
                                          * Initialize the local data
                                          * in the
@@ -492,6 +507,10 @@ namespace MeshWorker
                                          * assembled.
                                          */
        SmartPointer<MATRIX,MatrixSimple<MATRIX> > matrix;
+      /**
+       * A pointer to the object containing constraints.
+       */
+        SmartPointer<const ConstraintMatrix,MatrixSimple<MATRIX> > constraints;
        
                                         /**
                                          * The smallest positive
@@ -753,12 +772,13 @@ namespace MeshWorker
       /**
        * A pointer to the object containing the block structure.
        */
-      SmartPointer<const BlockInfo> block_info;
+      SmartPointer<const BlockInfo, 
+        MatrixLocalBlocksToGlobalBlocks<MATRIX, number> > block_info;
       /**
        * A pointer to the object containing constraints.
        */
       SmartPointer<const ConstraintMatrix, 
-       MatrixBlockVector<MATRIX> > constraints;
+       MatrixLocalBlocksToGlobalBlocks<MATRIX,number> > constraints;
       
                                         /**
                                          * The smallest positive
@@ -911,7 +931,7 @@ namespace MeshWorker
       /**
        * A pointer to the object containing the block structure.
        */
-      SmartPointer<const BlockInfo> block_info;
+      SmartPointer<const BlockInfo, MGMatrixLocalBlocksToGlobalBlocks<MATRIX, number> > block_info;
       
                                         /**
                                          * The smallest positive
@@ -1114,6 +1134,13 @@ namespace MeshWorker
       residuals = results;
     }
 
+    template <class VECTOR>
+    inline void
+    ResidualSimple<VECTOR>::initialize(const ConstraintMatrix& c)
+    {
+      constraints = &c;
+    }
+
     
     template <class VECTOR>
     template<int dim>
@@ -1130,8 +1157,16 @@ namespace MeshWorker
     ResidualSimple<VECTOR>::assemble(const DoFInfo<dim>& info)
     {
       for (unsigned int k=0;k<residuals.size();++k)
-       for (unsigned int i=0;i<info.vector(k).block(0).size();++i)
-         (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
+      {
+        if(constraints == 0)
+        {
+          for (unsigned int i=0;i<info.vector(k).block(0).size();++i)
+            (*residuals(k))(info.indices[i]) += info.vector(k).block(0)(i);
+        }
+        else
+          constraints->distribute_local_to_global(
+              info.vector(k).block(0), info.indices, (*residuals(k)));
+      }
     }
 
     
@@ -1143,10 +1178,20 @@ namespace MeshWorker
     {
       for (unsigned int k=0;k<residuals.size();++k)
        {
-         for (unsigned int i=0;i<info1.vector(k).block(0).size();++i)
-           (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
-         for (unsigned int i=0;i<info2.vector(k).block(0).size();++i)
-           (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
+          if(constraints == 0)
+          {
+            for (unsigned int i=0;i<info1.vector(k).block(0).size();++i)
+              (*residuals(k))(info1.indices[i]) += info1.vector(k).block(0)(i);
+            for (unsigned int i=0;i<info2.vector(k).block(0).size();++i)
+              (*residuals(k))(info2.indices[i]) += info2.vector(k).block(0)(i);
+          }
+          else
+          {
+            constraints->distribute_local_to_global(
+              info1.vector(k).block(0), info1.indices, (*residuals(k)));
+            constraints->distribute_local_to_global(
+              info2.vector(k).block(0), info2.indices, (*residuals(k)));
+          }
        }
     }
 
@@ -1189,20 +1234,20 @@ namespace MeshWorker
     {
         if(constraints == 0)
         {
-      for (unsigned int b=0;b<local.n_blocks();++b)
-       for (unsigned int j=0;j<local.block(b).size();++j)
-         {
-                                            // The coordinates of
-                                            // the current entry in
-                                            // DoFHandler
-                                            // numbering, which
-                                            // differs from the
-                                            // block-wise local
-                                            // numbering we use in
-                                            // our local vectors
-           const unsigned int jcell = this->block_info->local().local_to_global(b, j);
-           global(dof[jcell]) += local.block(b)(j);
-         }
+          for (unsigned int b=0;b<local.n_blocks();++b)
+            for (unsigned int j=0;j<local.block(b).size();++j)
+            {
+              // The coordinates of
+              // the current entry in
+              // DoFHandler
+              // numbering, which
+              // differs from the
+              // block-wise local
+              // numbering we use in
+              // our local vectors
+              const unsigned int jcell = this->block_info->local().local_to_global(b, j);
+              global(dof[jcell]) += local.block(b)(j);
+            }
         }
         else
           constraints->distribute_local_to_global(local, dof, global);
@@ -1251,6 +1296,14 @@ namespace MeshWorker
     {
       matrix = &m;
     }
+
+
+    template <class MATRIX>
+    inline void
+    MatrixSimple<MATRIX>::initialize(const ConstraintMatrix& c)
+    {
+      constraints = &c;
+    }
     
 
     template <class MATRIX >
@@ -1272,11 +1325,17 @@ namespace MeshWorker
     {
       AssertDimension(M.m(), i1.size());
       AssertDimension(M.n(), i2.size());
-      
+     
+     if(constraints == 0)
+     { 
       for (unsigned int j=0; j<i1.size(); ++j)
        for (unsigned int k=0; k<i2.size(); ++k)
          if (std::fabs(M(j,k)) >= threshold)
            matrix->add(i1[j], i2[k], M(j,k));
+     }
+     else
+       constraints->distribute_local_to_global(
+           M, i1, i2, *matrix);
     }
     
     
index 41f709fd65c0103a6625fd814f70c283fa72cc98..2a01c5839e4fe7d59e4c15e073fda0f28513c76f 100644 (file)
@@ -68,6 +68,18 @@ namespace MeshWorker
               bool gradients = false,
               bool hessians = false);
 
+                                      /**
+                                       * Does the same as the function
+                                        * above but it is possible to 
+                                        * select a block of the global
+                                        * vector.
+                                       */
+//      void add(const std::string& name,
+//               const unsigned int selected_block,
+//            bool values = true,
+//            bool gradients = false,
+//            bool hessians = false);
+
                                       /**
                                        * Initialize the selection
                                        * field with a data
@@ -449,6 +461,17 @@ namespace MeshWorker
   }
 
 
+  //inline void
+  //VectorSelector::add(const std::string& name,
+  //   const unsigned int block,
+  //   bool values, bool gradients, bool hessians)
+  //{
+  //  if (values) value_selection.add(name, block);
+  //  if (gradients) gradient_selection.add(name, block);
+  //  if (hessians) hessian_selection.add(name, block);  
+  //}
+
+
   template <typename DATA>
   inline void
   VectorSelector::initialize(const NamedData<DATA>& src)
index 5192221d94505856cd9016f6af18d60130f2b686..5d03cea5a495a504e5e5ee3a29d025d0cc0adb53 100644 (file)
@@ -82,6 +82,12 @@ namespace Algorithms
                                        * Read the parameters.
                                        */
       void initialize (ParameterHandler& param);
+
+                                      /**
+                                       * Initialize the pointer 
+                                        * data_out for debugging.
+                                       */
+      void initialize (OutputOperator<VECTOR>& output);
       
                                       /**
                                        * The actual Newton
@@ -116,13 +122,20 @@ namespace Algorithms
                                       /**
                                        * The operator computing the residual.
                                        */
-      SmartPointer<Operator<VECTOR> > residual;
+      SmartPointer<Operator<VECTOR>, Newton<VECTOR> > residual;
 
                                       /**
                                        * The operator applying the
                                        * inverse derivative to the residual.
                                        */
-      SmartPointer<Operator<VECTOR> > inverse_derivative;
+      SmartPointer<Operator<VECTOR>, Newton<VECTOR> > inverse_derivative;
+
+                                      /**
+                                       * The operator handling the output 
+                                        * in case the debug_vectors is true.
+                                        * Call the initialize function first.
+                                       */
+      SmartPointer<OutputOperator<VECTOR>, Newton<VECTOR> > data_out;
       
                                       /**
                                        * This flag is set by the
@@ -166,6 +179,8 @@ namespace Algorithms
                                        * every Newton step.
                                        */
       double assemble_threshold;
+
+    public:
                                       /**
                                        * Print residual, update and
                                        * updated solution after each
index 6f1e0173f92e006034e4e6b911764bd743a9cb82..3728c1fe8c7b49b48e17789757e3ab32b2b803d6 100644 (file)
@@ -60,6 +60,12 @@ namespace Algorithms
     param.leave_subsection ();
   }
 
+  template <class VECTOR>
+  void
+  Newton<VECTOR>::initialize (OutputOperator<VECTOR>& output)
+  {
+    data_out = &output;  
+  }
 
   template <class VECTOR>
   void
@@ -125,6 +131,16 @@ namespace Algorithms
                    << e.last_residual << std::endl;
          }
 
+       if (debug_vectors)
+         {
+           NamedData<VECTOR*> out;
+            VECTOR* p = &u; out.add(p, "solution");
+           p = Du; out.add(p, "update");
+           p = res; out.add(p, "residual");
+            *data_out << step;
+           *data_out << out;
+         }
+
        u.add(-1., *Du);
        old_residual = resnorm;
        (*residual)(out1, src1);
@@ -147,19 +163,6 @@ namespace Algorithms
            (*residual)(out1, src1);
            resnorm = res->l2_norm();
          }
-
-       if (debug_vectors)
-         {
-           std::ostringstream streamOut;
-           streamOut << "Newton_" << std::setw(3) << std::setfill('0') << step;
-           std::string name = streamOut.str();
-           NamedData<VECTOR*> out;
-           out.add(&u, "solution");
-           out.add(Du, "update");
-           out.add(res, "residual");
-//TODO: Implement!
-//         app->data_out(name, out);
-         }
       }
     deallog.pop();
 
index fd6d6ae98034c25411f8d3d9b8021c5a874c0721..cad77faf3b5e1f375486e0d8cce2ee30e0b3a7de 100644 (file)
@@ -16,6 +16,9 @@
 #include <base/config.h>
 #include <base/named_data.h>
 #include <base/event.h>
+#include <dofs/dof_handler.h>
+
+#include <fstream>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -94,11 +97,21 @@ namespace Algorithms
   template <class VECTOR>
   class OutputOperator : public Subscriptor
   {
+    OutputOperator(const OutputOperator<VECTOR>&);
     public:
+    OutputOperator ();
                                       /**
                                        * Empty virtual destructor.
                                        */
       virtual ~OutputOperator();
+      /**
+       * Set the stream #os to
+       * which data is written. If
+       * no stream is selected with
+       * this function, data goes
+       * to #deallog.
+       */
+      void initialize_stream(std::ostream& stream);
                                       /**
                                        * Set the current step.
                                        */
@@ -107,9 +120,11 @@ namespace Algorithms
                                       /**
                                        * Output all the vectors in NamedData.
                                        */
-      virtual OutputOperator<VECTOR>& operator<< (const NamedData<VECTOR*> vectors);
+      virtual OutputOperator<VECTOR>& operator<< (const NamedData<VECTOR*>& vectors);
     protected:
       unsigned int step;
+    private:
+      std::ostream* os;
   };
 
   template <class VECTOR>
@@ -119,6 +134,27 @@ namespace Algorithms
     step = s;
     return *this;
   }
+
+  template <class VECTOR, int dim, int spacedim=dim>
+  class DoFOutputOperator : public OutputOperator<VECTOR>
+  {
+    public:
+      void initialize (DoFHandler<dim, spacedim>& dof_handler);
+
+      virtual OutputOperator<VECTOR>& 
+        operator<<(const NamedData<VECTOR*>& vectors);
+
+    private:
+      SmartPointer<DoFHandler<dim, spacedim>,
+        DoFOutputOperator<VECTOR, dim, spacedim> > dof;
+  };
+
+  template <class VECTOR, int dim, int spacedim>
+  void DoFOutputOperator<VECTOR, dim, spacedim>::initialize(
+      DoFHandler<dim, spacedim>& dof_handler)
+  {
+    dof = &dof_handler;  
+  }
 }
 
 
index 20b7bc96dae5460b5f18ec5bfdc098ff021917af..a2ac9e57f633279129a49e62cf12f73c70d5752d 100644 (file)
@@ -11,6 +11,7 @@
 //---------------------------------------------------------------------------
 
 #include <numerics/operator.h>
+#include <numerics/data_out.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -19,6 +20,7 @@ namespace Algorithms
   template <class VECTOR>
   Operator<VECTOR>::~Operator()
   {}
+
   
   
   template <class VECTOR>
@@ -40,20 +42,58 @@ namespace Algorithms
   OutputOperator<VECTOR>::~OutputOperator()
   {}
   
+  template <class VECTOR>
+  OutputOperator<VECTOR>::OutputOperator()
+  :
+    os(0)
+  {}
+
+  template <class VECTOR>
+  void OutputOperator<VECTOR>::initialize_stream(std::ostream& stream)
+  {
+    os =&stream;
+  }
   
   template <class VECTOR>
   OutputOperator<VECTOR>&
-  OutputOperator<VECTOR>::operator<< (const NamedData<VECTOR*> vectors)
+  OutputOperator<VECTOR>::operator<< (const NamedData<VECTOR*>& vectors)
   {
-    std::cout << ' ' << step;
-    for (unsigned int i=0;i<vectors.size();++i)
-      vectors(i)->print(std::cout);
-    std::cout << std::endl;
+    if(os == 0)
+    {
+      //TODO: make this possible
+      //deallog << ' ' << step;
+      //for (unsigned int i=0;i<vectors.size();++i)
+      //  vectors(i)->print(deallog);
+      //deallog << std::endl;
+    }
+    else
+    {
+      (*os) << ' ' << step;
+      for (unsigned int i=0;i<vectors.size();++i)
+        vectors(i)->print(*os);
+      (*os) << std::endl;
+    }
+    return *this;
+  }
+
+  template <class VECTOR, int dim, int spacedim>
+  OutputOperator<VECTOR>& 
+  DoFOutputOperator<VECTOR, dim, spacedim>::operator<<(
+      const NamedData<VECTOR*>& vectors)
+  {
+    Assert ((dof!=0), ExcNotInitialized());
+    DataOut<dim> out;
+    out.attach_dof_handler (*dof);
+    out.add_data_vector (*vectors(vectors.find("solution")), "solution");
+    out.add_data_vector (*vectors(vectors.find("update")), "update");
+    out.add_data_vector (*vectors(vectors.find("residual")), "residual");
+    std::ostringstream streamOut;
+    streamOut << "Newton_" << std::setw(3) << std::setfill('0') << this->step;
+    std::ofstream out_filename (streamOut.str().c_str());
+    out.build_patches (2);
+    out.write_gnuplot (out_filename);
     return *this;
   }
-  
-  
-  
 }
 
 DEAL_II_NAMESPACE_CLOSE
index f23595084ef4bdf6a6e177fac94511b29f4324b9..668eeb2a6d8d0b0570e884185f994dd90d734413 100644 (file)
@@ -1130,7 +1130,9 @@ add_data_vector<V > (const V                  &, \
     INSTANTIATE(DH,D0,D1,D2,D3,Vector<double>); \
     INSTANTIATE(DH,D0,D1,D2,D3,Vector<float>); \
     INSTANTIATE(DH,D0,D1,D2,D3,BlockVector<double>) ; \
-    INSTANTIATE(DH,D0,D1,D2,D3,BlockVector<float>)
+    INSTANTIATE(DH,D0,D1,D2,D3,BlockVector<float>); \
+    INSTANTIATE(DH,D0,D1,D2,D3,Vector<long double>); \
+    INSTANTIATE(DH,D0,D1,D2,D3,BlockVector<long double>)
 #else
 # define INSTANTIATE_VECTORS(DH,D0,D1,D2,D3)    \
     INSTANTIATE(DH,D0,D1,D2,D3,Vector<double>); \
index fe946a62567f2e54acfcbade8512a9e0eb33f27b..90a137b58fe320a85e404ef303e42621b6d0e086 100644 (file)
@@ -31,6 +31,6 @@ DEAL_II_NAMESPACE_OPEN
 
 using namespace Algorithms;
 
-#include "operator.inst"
+#include "operator.all_dimensions.inst"
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/deal.II/source/numerics/operator.all_dimensions.inst.in b/deal.II/deal.II/source/numerics/operator.all_dimensions.inst.in
new file mode 100644 (file)
index 0000000..d4cf0f7
--- /dev/null
@@ -0,0 +1,20 @@
+//---------------------------------------------------------------------------
+//    $Id: vectors.inst.in 19046 2009-07-08 19:30:23Z bangerth $
+//    Version: $Name$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+for (VEC : SERIAL_VECTORS)
+{  
+  template class Operator<VEC>;
+  template class OutputOperator<VEC>;
+  template class Newton<VEC>;
+  template class ThetaTimestepping<VEC>;
+}
diff --git a/deal.II/deal.II/source/numerics/operator.cc b/deal.II/deal.II/source/numerics/operator.cc
new file mode 100644 (file)
index 0000000..fe946a6
--- /dev/null
@@ -0,0 +1,36 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <base/parameter_handler.h>
+#include <base/logstream.h>
+#include <lac/vector_memory.h>
+
+#include <numerics/operator.templates.h>
+#include <numerics/newton.templates.h>
+#include <numerics/theta_timestepping.templates.h>
+
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/petsc_vector.h>
+#include <lac/petsc_block_vector.h>
+#include <lac/trilinos_vector.h>
+#include <lac/trilinos_block_vector.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+using namespace Algorithms;
+
+#include "operator.inst"
+
+DEAL_II_NAMESPACE_CLOSE
index 7b5086980fd9d74323c5fbe6e94931f44b11c79f..631ec8e98387e6bb9c60000feffacbdda4d7e2ed 100644 (file)
@@ -13,7 +13,5 @@
 
 for (VEC : SERIAL_VECTORS)
 {  
-  template class Operator<VEC>;
-  template class Newton<VEC>;
-  template class ThetaTimestepping<VEC>;
+  template class DoFOutputOperator<VEC, deal_II_dimension, deal_II_dimension>;
 }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.