oomph-lib

00001 //LIC// ====================================================================
00002 //LIC// This file forms part of oomph-lib, the object-oriented, 
00003 //LIC// multi-physics finite-element library, available 
00004 //LIC// at http://www.oomph-lib.org.
00005 //LIC// 
00006 //LIC//           Version 0.90. August 3, 2009.
00007 //LIC// 
00008 //LIC// Copyright (C) 2006-2009 Matthias Heil and Andrew Hazel
00009 //LIC// 
00010 //LIC// This library is free software; you can redistribute it and/or
00011 //LIC// modify it under the terms of the GNU Lesser General Public
00012 //LIC// License as published by the Free Software Foundation; either
00013 //LIC// version 2.1 of the License, or (at your option) any later version.
00014 //LIC// 
00015 //LIC// This library is distributed in the hope that it will be useful,
00016 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 //LIC// Lesser General Public License for more details.
00019 //LIC// 
00020 //LIC// You should have received a copy of the GNU Lesser General Public
00021 //LIC// License along with this library; if not, write to the Free Software
00022 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
00023 //LIC// 02110-1301  USA.
00024 //LIC// 
00025 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
00026 //LIC// 
00027 //LIC//====================================================================
00028 #include "orthpoly.h"
00029 
00030 namespace oomph
00031 {
00032 
00033 namespace Orthpoly
00034 {
00035 // Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1
00036 void gll_nodes(const unsigned &Nnode, Vector<double> &x)
00037 {
00038   double z,zold,del;
00039   unsigned p=Nnode-1;
00040   x.resize(Nnode);
00041   if(Nnode<2) 
00042    {
00043     throw OomphLibError("Invalid number of nodes",
00044                         "Orthpoly::gll_nodes()",
00045                         OOMPH_EXCEPTION_LOCATION);
00046    }
00047   else if (Nnode==2){
00048     x[0]=-1.0;x[1]=1.0;
00049   }
00050   else if (Nnode==3){
00051     x[0]=-1;x[1]=0.0;x[2]=1.0;
00052   }
00053   else if (Nnode==4){
00054     x[0]=-1.0;x[1]=-std::sqrt(5.0)/5.0;x[2]=-x[1];x[3]=1.0;
00055   }
00056   else
00057     {
00058       unsigned mid;
00059       if(p%2==0){
00060         mid = p/2;
00061         x[mid]=0.0;
00062       }
00063       else mid=p/2+1;
00064       for(unsigned j=1;j<mid;j++) 
00065         {
00066           z=-std::cos(j*pi/double(p));
00067           do{
00068            del= dlegendre(p,z)/ddlegendre(p,z);
00069            zold=z;
00070            z=zold-del;
00071           } while(std::abs(z-zold) > eps);
00072           x[j]=z;
00073           x[p-j]=-z;
00074         }
00075       x[0]=-1.0;x[p]=1.0;
00076     }
00077 }
00078 
00079 // This version of gll_nodes calculates the abscissas AND weights
00080 void gll_nodes(const unsigned &Nnode, Vector<double> &x, Vector<double> &w)
00081 {
00082  gll_nodes(Nnode,x);
00083  // Now calculate the corresponding weights
00084  double l_z;
00085  unsigned p = Nnode-1;
00086  for (unsigned i=0;i<p+1;i++){
00087   l_z=legendre(p,x[i]);
00088   w[i]=2.0/(p*(p+1)*l_z*l_z);
00089  }
00090 }
00091 
00092 
00093 // Calculates the Gauss Legendre abscissas of degree p=Nnode-1
00094 void gl_nodes(const unsigned &Nnode, Vector<double> &x)
00095 {
00096  double z,zold,del;
00097  unsigned p=Nnode-1;
00098  x.resize(Nnode);
00099  if(Nnode<2)
00100   {
00101    throw OomphLibError("Invalid number of nodes",
00102                        "Orthpoly::gl_nodes()",
00103                        OOMPH_EXCEPTION_LOCATION);
00104   }
00105  else if (Nnode==2){
00106   x[0]=-1.0/3.0*std::sqrt(3.0);x[1]=-x[0];
00107  }
00108  else
00109   {
00110    unsigned mid;
00111    if(p%2==0){
00112     mid = p/2;
00113     x[mid]=0.0;
00114    }
00115    else mid=p/2+1;
00116    for(unsigned j=0;j<mid;j++) 
00117     {
00118      z=-std::cos((2.0*j+1.0)*pi/(2*double(p)+2.0));
00119      do{
00120       del = legendre(p+1,z)/dlegendre(p+1,z);
00121       zold=z;
00122             z=zold-del;
00123      } while(std::abs(z-zold)>eps);
00124      x[j]=z;
00125      x[p-j]=-z;
00126     }
00127   }
00128 }
00129 
00130 // This version of gl_nodes calculates the abscissas AND weights
00131 void gl_nodes(const unsigned &Nnode, Vector<double> &x, Vector<double> &w)
00132 {
00133  gl_nodes(Nnode,x);
00134   // Now calculate the corresponding weights
00135   double dl_z;
00136   unsigned p = Nnode-1;
00137   for (unsigned i=0;i<p+1;i++){
00138    dl_z=dlegendre(p+1,x[i]);
00139    w[i]=2.0/(1-x[i]*x[i])/(dl_z*dl_z);
00140   }
00141 }
00142 }
00143 
00144 }
orthpoly.cc
