matrix is long double

git-svn-id: https://svn.dealii.org/trunk@5509 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/contrib/utilities/Makefile.in b/deal.II/contrib/utilities/Makefile.in
index 72c1a971415c0dc60479dfda53db24aa3d46f52f..2a7bf3a21d3dbeb31a2bbeaeaab784188c89fd46 100644 (file)
--- a/deal.II/contrib/utilities/Makefile.in
+++ b/deal.II/contrib/utilities/Makefile.in
@@ -1,6 +1,6 @@
 ############################################################
 # $Id$
-# Copyright (C) 2000, 2001 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002 by the deal.II authors
 ############################################################
 
 ############################################################
@@ -9,31 +9,31 @@
 
 D = @DEAL2_DIR@
 include $D/common/Make.global_options
-debug-mode = on
+debug-mode = off
 
-libraries = $(lib-deal2-1d.g) \
-            $(lib-deal2-2d.g) \
-            $(lib-deal2-3d.g) \
-            $(lib-lac.g)  \
-            $(lib-base.g)
+libraries = $(lib-deal2-1d.o) \
+            $(lib-deal2-2d.o) \
+            $(lib-deal2-3d.o) \
+            $(lib-lac.o)  \
+            $(lib-base.o)
 
 ############################################################
 # First how to create executables, including all necessary
 # flags:
 ############################################################
 
-flags     = $(CXXFLAGS.g)
+flags     = $(CXXFLAGS.o)
 
 ifeq ($(findstring gcc,$(GXX_VERSION)),gcc)
 flags += -Wno-missing-noreturn
 endif
 
-%.o : %.cc Makefile
-       @echo =====optimized===== $<
-       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
 %.g.o : %.cc Makefile
        @echo =====debug========= $<
        @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+%.o : %.cc Makefile
+       @echo =====optimized===== $<
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
 %.exe :
        @echo =====linking======= $@
        @$(CXX) $(LDFLAGS) -o $@ $^ $(LIBS)

diff --git a/deal.II/contrib/utilities/embedding.cc b/deal.II/contrib/utilities/embedding.cc
index f04b20bf23fe9347b50875b0af46dd7f851b851f..1cbc34ed101e72debbaeace1129e46f928677f3e 100644 (file)
--- a/deal.II/contrib/utilities/embedding.cc
+++ b/deal.II/contrib/utilities/embedding.cc
@@ -14,8 +14,9 @@
 #include <grid/grid_generator.h>
 #include <fe/fe_q.h>
 #include <fe/fe_dgq.h>
+#include <fe/fe_dgp.h>
 #include <fe/fe_system.h>
-#include <fe/mapping_q1.h>
+#include <fe/mapping_cartesian.h>
 #include <fe/fe_values.h>
 #include <vector>
 #include <iomanip>
@@ -23,22 +24,32 @@
 #include <strstream>
 #include <string>
 
+/*
+ * Enter name of the finite element family here.
+ */
+
+#define ELNAME "FE_DGP"
+#define PUSH(k) FE_DGP<dim> q ## k (k);\
+ elements.push_back (ElPair(&q ## k, "dgp" # k))
+
+
 
 template <int dim>
-void compute_embedding (Triangulation<dim>& tr_coarse,
+void compute_embedding (unsigned int degree,
+                       Triangulation<dim>& tr_coarse,
                        Triangulation<dim>& tr_fine,
                        const FiniteElement<dim>& fe_coarse,
                        const FiniteElement<dim>& fe_fine,
                        const char* name)
 {
-  cerr << name << '<' << dim << '>' << endl;
+  deallog.push(name);
   
   DoFHandler<dim> dof_coarse(tr_coarse);
   dof_coarse.distribute_dofs(fe_coarse);
   DoFHandler<dim> dof_fine(tr_fine);
   dof_fine.distribute_dofs(fe_fine);
   
-  FullMatrix<double> A(dof_fine.n_dofs());
+  FullMatrix<long double> A(dof_fine.n_dofs());
   Vector<long double> f(dof_fine.n_dofs());
   Vector<long double> u(dof_fine.n_dofs());
   vector<vector<vector<double> > >
@@ -50,8 +61,8 @@ void compute_embedding (Triangulation<dim>& tr_coarse,
     = dof_coarse.begin_active();
   vector<unsigned int> indices(fe_fine.dofs_per_cell);
 
-  MappingQ1<dim> mapping;
-  QGauss<dim> q_fine(12);
+  MappingCartesian<dim> mapping;
+  QGauss<dim> q_fine(degree+1);
   
   FEValues<dim> fine (mapping, fe_fine, q_fine,
                      update_q_points
@@ -98,7 +109,7 @@ void compute_embedding (Triangulation<dim>& tr_coarse,
 //      A.print_formatted(cout, 2, false, 4, "~", 9);
       FullMatrix<double> P(A.n());
       P.invert(A);
-      SolverControl control (100, 1.e-20, true, false);
+      SolverControl control (100, 1.e-24, true, false);
       PrimitiveVectorMemory<Vector<long double> > mem;
       SolverRichardson<Vector<long double> > solver(control, mem);
       solver.solve(A,u,f,P);
@@ -117,10 +128,9 @@ void compute_embedding (Triangulation<dim>& tr_coarse,
 
   for (unsigned int cell=0;cell<tr_fine.n_active_cells();++cell)
     {
-      cout << "static double[] "
+      cout << "static const double "
           << name << "_"
-          << dim << "d_"
-          << cell << " =" << endl << '{' << endl;
+          << cell << "[] =" << endl << '{' << endl;
       for (unsigned int i=0;i<fe_fine.dofs_per_cell;++i)
        {
          for (unsigned int j=0;j<fe_coarse.dofs_per_cell;++j)
@@ -131,17 +141,19 @@ void compute_embedding (Triangulation<dim>& tr_coarse,
        }
       cout << "};" << endl << endl;
     }
+  deallog.pop();
 }
 
-#define PUSH_Q(k) FE_Q<dim> q ## k (k);\
- elements.push_back (ElPair(&q ## k, "q" # k))
-
-#define PUSH_DGQ(k) FE_DGQ<dim> dgq ## k(k);\
- elements.push_back (ElPair(&dgq ## k, "dgq" # k))
 
 template <int dim>
 void loop ()
 {
+  char prefix[3];
+  sprintf(prefix, "%dd", dim);
+  deallog.push(prefix);
+  
+  cout << "namespace " << ELNAME << "_" << dim << "d\n{";
+
   Triangulation<dim> tr_coarse;
   Triangulation<dim> tr_fine;
   GridGenerator::hyper_cube (tr_coarse, 0, 1);
@@ -151,37 +163,43 @@ void loop ()
   typedef pair<const FiniteElement<dim>*, const char*> ElPair ;
   vector <ElPair> elements(0);
 
-        PUSH_DGQ(0);
-        PUSH_DGQ(1);
-        PUSH_DGQ(2);
-//        PUSH_DGQ(3);
-        PUSH_DGQ(5);
-        PUSH_DGQ(6);
-        PUSH_DGQ(7);
-      //  PUSH_Q(0);
-      PUSH_Q(1);
-      PUSH_Q(2);
-      PUSH_Q(3);
-      PUSH_Q(4);
+                                  /*
+                                   * List element degrees for
+                                   * computation here.
+                                   */
+  PUSH(0);
+  PUSH(1);
+  PUSH(2);
+  PUSH(3);
+  PUSH(4);
+  PUSH(5);
+  PUSH(6);
 
   char* name = new char[100];
 
-                                  // Embed all lower spaces into higher
+                                  // Embed all lower spaces into
+                                  // higher or just the same degree
+                                  // on different grids.
+  bool same_degree_only = true;
+  
   unsigned int n = elements.size();
   for (unsigned int i=0;i<n;++i)
-    for (unsigned int j=i;j<=i;++j)
+    for (unsigned int j=((same_degree_only) ? i : 0);j<=i;++j)
       {
        ostrstream os (name, 99);
        os << elements[i].second << "_into_"
           << elements[j].second << "_refined" << ends;
        
-       compute_embedding (tr_coarse, tr_fine,
+       compute_embedding (i, tr_coarse, tr_fine,
                           *elements[i].first,
                           *elements[j].first,
                           name);
       }
 
   delete [] name;
+  cout << "};\n";
+  deallog.pop();
+  
 }
 
 int main ()

diff --git a/deal.II/contrib/utilities/simplify.pl b/deal.II/contrib/utilities/simplify.pl
new file mode 100644 (file)
index 0000000..f6f97b1
--- /dev/null
+++ b/deal.II/contrib/utilities/simplify.pl
@@ -0,0 +1,8 @@
+s! -?\d\.\d+e-1[0123456789]/27.! 0.!g;
+s! 0/27.! 0.!g;
+s! (-?)27/27.! ${1}1.!g;
+s! (-?)13\.5/27.! ${1}.5!g;
+s! (-?)20\.25/27.! ${1}.75!g;
+s! (-?)6\.75/27.! ${1}.25!g;
+s! (-?)10\.125/27.! ${1}.375!g;
+s! (-?)3\.375/27.! ${1}.125!g;