]> https://gitweb.dealii.org/ - dealii.git/commitdiff
First implementation of values and grads. grad_grads not yet implemented.
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Mon, 18 Dec 2000 10:48:21 +0000 (10:48 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Mon, 18 Dec 2000 10:48:21 +0000 (10:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@3548 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/tensor_product_polynomials.cc [new file with mode: 0644]

diff --git a/deal.II/base/source/tensor_product_polynomials.cc b/deal.II/base/source/tensor_product_polynomials.cc
new file mode 100644 (file)
index 0000000..c7d7553
--- /dev/null
@@ -0,0 +1,266 @@
+//----------------------  tensor_product_polynomials.cc  ------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000 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.
+//
+//----------------------  tensor_product_polynomials.cc  ------------
+
+
+#include <base/exceptions.h>
+#include <base/tensor_product_polynomials.h>
+
+#include <math.h>
+
+
+template <int dim>
+TensorProductPolynomials<dim>::TensorProductPolynomials(
+  const vector<SmartPointer<Polynomial> > &pols):
+               polynomials(pols)
+{}
+
+                                // first version of this function
+//  template <int dim>
+//  void TensorProductPolynomials<dim>::shape_values_and_grads(
+//    const Point<dim> &p,
+//    vector<double> &values,
+//    vector<Tensor<1,dim> > &grads,
+//    vector<Tensor<2,dim> > &grad_grads) const
+//  {
+//    unsigned int n_pols=polynomials->size();
+//    unsigned int n_pols2=n_pols*n_pols;
+//    unsigned int n_tensor_pols=pow(n_pols, dim);
+//    Assert(values.size()==n_tensor_pols || values.size()==0,
+//      ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
+//    Assert(grads.size()==n_tensor_pols|| grads.size()==0,
+//      ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
+//    Assert(grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
+//      ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
+
+//    unsigned int v_size=0;
+//    bool update_values=false, update_grads=false, update_grad_grads=false;
+//    if (values.size()==n_tensor_pols)
+//      {
+//        update_values=true;
+//        v_size=1;
+//      }
+//    if (grads.size()==n_tensor_pols)
+//      {
+//        update_grads=true;
+//        v_size=2;
+//      }
+//    if (grad_grads.size()==n_tensor_pols)
+//      {
+//        update_grad_grads=true;
+//        v_size=3;
+//      }
+
+//    vector<vector<vector<double> > > v(
+//      dim, vector<vector<double> > (n_pols, vector<double> (v_size)));
+
+//    for (unsigned int d=0; d<dim; ++d)
+//      {
+//        vector<vector<double> > &v_d=v[d];
+//        Assert(v_d.size()==n_pols, ExcInternalError());
+//        for (unsigned int i=0; i<n_pols; ++i)
+//     polynomials->operator[](i)->value(p(d), v_d[i]);
+//      }
+  
+//    if (update_values)
+//      {
+//                                    // values[i]=p_j(x)
+//        for (unsigned i=0; i<n_tensor_pols;)
+//     for (unsigned int j=0; j<n_pols; ++j, ++i)
+//       values[i]=v[0][j][0];
+
+//                                    // values[i]*=p_j(y)
+//        if (dim>1)
+//          for (unsigned i=0; i<n_tensor_pols;)
+//       for (unsigned j=0; j<n_pols; ++j)
+//         {
+//           double pjy=v[1][j][0];
+//           for (unsigned int k=0; k<n_pols; ++k, ++i)
+//             values[i]*=pjy;
+//         }
+      
+//                                    // values[i]*=p_j(z)
+//        if (dim>2)
+//     for (unsigned int i=0; i<n_tensor_pols;)
+//       for (unsigned int j=0; j<n_pols; ++j)
+//         {
+//           double pjz=v[2][j][0];
+//           for (unsigned int k=0; k<n_pols2; ++k, ++i)
+//             values[i]*=pjz;
+//         }
+      
+//        Assert(dim<=3, ExcNotImplemented());
+//      }
+
+      
+//    if (update_grads)
+//      {
+//                                    // grads[i][0]=dp_j(x)
+//                                    // grads[i][1]=p_j(x)
+//                                    // grads[i][2]=p_j(x)
+//        for (unsigned i=0; i<n_tensor_pols;)
+//     for (unsigned int j=0; j<n_pols; ++j, ++i)
+//       switch (dim)
+//         {
+//           case 3: grads[i][2]=v[0][j][0];
+//           case 2: grads[i][1]=v[0][j][0];
+//           case 1: grads[i][0]=v[0][j][1]; break;
+//           default: Assert(false, ExcNotImplemented());
+//         }
+
+
+//                                    // grad[i][0]*=p_j(y)
+//                                    // grad[i][1]*=dp_j(y)
+//                                    // grad[i][2]*=p_j(y)
+//        if (dim>1)
+//     for (unsigned i=0; i<n_tensor_pols;)
+//       for (unsigned j=0; j<n_pols; ++j)
+//         {
+//           double pjy=v[1][j][0];
+//           double dpjy=v[1][j][1];
+//           for (unsigned int k=0; k<n_pols; ++k, ++i)
+//             {
+//               grads[i][0]*=pjy;
+//               grads[i][1]*=dpjy;
+//               if (dim==3) grads[i][2]*=pjy;
+//               Assert(dim<=3, ExcNotImplemented());
+//             }
+//         }
+      
+//                                    // grad[i][0]*=p_j(z)
+//                                    // grad[i][1]*=p_j(z)
+//                                    // grad[i][2]*=dp_j(z)
+//        if (dim>2)
+//     for (unsigned i=0; i<n_tensor_pols;)
+//       for (unsigned j=0; j<n_pols; ++j)
+//         {
+//           double pjz=v[2][j][0];
+//           double dpjz=v[2][j][1];
+//           for (unsigned int k=0; k<n_pols2; ++k, ++i)
+//             {
+//               grads[i][0]*=pjz;
+//               grads[i][1]*=pjz;
+//               grads[i][2]*=dpjz;
+//               Assert(dim<=3, ExcNotImplemented());
+//             }
+//         }
+//      }
+
+//    if (update_grad_grads)
+//      {
+//        Assert(false, ExcNotImplemented());
+//      }
+  
+//  }
+
+
+
+template <typename number>
+number pow(number x, unsigned int y)
+{
+  number value=1;
+  for (unsigned int i=0; i<y; ++i)
+    value*=x;
+  return value;
+}
+
+
+                                // second and shorter version
+template <int dim>
+void TensorProductPolynomials<dim>::shape_values_and_grads(
+  const Point<dim> &p,
+  vector<double> &values,
+  vector<Tensor<1,dim> > &grads,
+  vector<Tensor<2,dim> > &grad_grads) const
+{
+  unsigned int n_pols=polynomials.size();
+  unsigned int n_tensor_pols=pow(n_pols, dim);
+  Assert(values.size()==n_tensor_pols || values.size()==0,
+        ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
+  Assert(grads.size()==n_tensor_pols|| grads.size()==0,
+        ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
+  Assert(grad_grads.size()==n_tensor_pols|| grad_grads.size()==0,
+        ExcDimensionMismatch2(grad_grads.size(), n_tensor_pols, 0));
+
+  unsigned int v_size=0;
+  bool update_values=false, update_grads=false, update_grad_grads=false;
+  if (values.size()==n_tensor_pols)
+    {
+      update_values=true;
+      v_size=1;
+    }
+  if (grads.size()==n_tensor_pols)
+    {
+      update_grads=true;
+      v_size=2;
+    }
+  if (grad_grads.size()==n_tensor_pols)
+    {
+      update_grad_grads=true;
+      v_size=3;
+    }
+
+  vector<vector<vector<double> > > v(
+    dim, vector<vector<double> > (n_pols, vector<double> (v_size)));
+
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      vector<vector<double> > &v_d=v[d];
+      Assert(v_d.size()==n_pols, ExcInternalError());
+      for (unsigned int i=0; i<n_pols; ++i)
+       polynomials[i]->value(p(d), v_d[i]);
+    }
+  
+  vector<unsigned int> n_pols_to(dim+1);
+  n_pols_to[0]=1;
+  for (unsigned int i=0; i<dim; ++i)
+    n_pols_to[i+1]=n_pols_to[i]*n_pols;
+  Assert(n_pols_to[dim]==n_tensor_pols, ExcInternalError());
+  
+  if (update_values)
+    {
+      for (unsigned int i=0; i<n_tensor_pols; ++i)
+       values[i]=1;
+      
+      for (unsigned int d=0; d<dim; ++d)
+       for (unsigned int i=0; i<n_tensor_pols; ++i)
+         values[i]*=v[d][(i/n_pols_to[d])%n_pols][0];
+    }
+
+      
+  if (update_grads)
+    {
+                                      // is Tensor<1,dim> &gradsi=grads[i]; grads[d];
+                                      // faster than grads[i][d]??
+      for (unsigned int i=0; i<n_tensor_pols; ++i)
+       for (unsigned int d=0; d<dim; ++d)
+         grads[i][d]=1.;
+
+      for (unsigned int x=0; x<dim; ++x)
+       for (unsigned int i=0; i<n_tensor_pols; ++i)
+         for (unsigned int d=0; d<dim; ++d)
+           grads[i][d]*=v[x][(i/n_pols_to[x])%n_pols][d==x];
+    }
+
+  if (update_grad_grads)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+  
+}
+
+
+
+  
+template class TensorProductPolynomials<1>;
+template class TensorProductPolynomials<2>;
+template class TensorProductPolynomials<3>;

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.