merged siac

library/CMakeLists.txt

 SET(LibrarySubDirs FieldUtils GlobalMapping LibUtilities LocalRegions
-    Collections MultiRegions SpatialDomains StdRegions)
+    Collections MultiRegions SpatialDomains StdRegions LSIAC LSIAC2D)
 SET(UnitTestSubDirs UnitTests)
 SET(DemoSubDirs     Demos)
 SET(TimingsSubDirs  Timings)

library/LSIAC/BSplines.h

+#ifndef LSIAC_BSPLINES_H
+#define LSIAC_BSPLINES_H
+
+#include <iostream>
+#include <vector>
+#include "NektarBaseClass.h"
+using namespace std;
+
+/// This class is the base class for all the B-Splines. 
+/** This class is useful when dynmaically creating object of its subclasses.
+	All the BSplines need by the fitler would be subclasses of this filter.
+*/ 
+
+namespace Nektar
+{
+namespace LSIAC
+{
+
+
+class BSplines: public NektarBaseClass{
+	protected:
+	public:
+};
+/*
+class GeneralBSplines: public BSplines{
+    // data
+	public:
+		vector<double>* m_knotVector;
+		int m_Order;
+
+	// functions
+	public:
+	GeneralBSplines(const vector<double> knots,const int order)
+	{
+		m_knotVector = new vector<double>(0);
+		SetKnotVector(knots);
+		SetOrder(order);
+	}
+
+	bool SetKnotVector(vector<double> knots)
+	{
+		this->m_knotVector->clear();
+		this->m_knotVector->resize(knots.size());
+		this->m_knotVector->assign(knots.begin(),knots.end());
+		return true;
+	}
+
+	bool GetKnotVector ( vector<double> &knots) 
+	{
+		knots = *m_knotVector;
+		return true; }
+	
+	bool SetOrder(int order)
+	{
+		m_Order = order;
+		return true; }
+
+	int  GetOrder() const
+	{return m_Order; }
+
+	bool EvaluateBSplines(const  vector<double> t_pos,const vector<double> knots,
+					vector<double> &t_values)
+	{return true; }
+
+	bool EvaluateBSplines (const vector<double> t_pos, vector<double> &t_values)const
+	{return true; }
+
+};
+
+int main()
+{
+ 	cout << " Enterend into main functions " << endl;
+	double mydoubles[] = { 0.0, 1.0, 2.0, 3.0, 4.0};	
+	vector<double> knots(mydoubles, mydoubles + sizeof(mydoubles)/sizeof(double));	
+	GeneralBSplines gsp(knots, 3);	
+	cout << gsp.m_knotVector->at(3)<< endl; 
+	cout <<  knots.size()<< endl;
+	vector<double> *knotsout;
+	gsp.GetKnotVector(*knotsout);
+	knotsout->at(3) = 10;
+	cout << knotsout->size()<< endl;
+	cout << knotsout->at(0)<< endl; 
+	cout << knotsout->at(1)<< endl; 
+	cout << knotsout->at(2)<< endl; 
+	cout << knotsout->at(3)<< endl; 
+	cout << gsp.m_knotVector->at(3)<< endl; 
+	
+	return 1;
+}
+*/
+
+}
+}
+
+
+#endif

library/LSIAC/CMakeLists.txt

+SET(LSIACHeaders
+	NektarBaseClass.h
+	BSplines.h
+	FilterType.h
+	GeneralBSplines.h
+	CentralBSplines.h
+	SIACFilter.h
+	SymmetricSIAC.h
+	OneSidedSIAC.h
+	HandleMesh.h
+	HandleNekMesh.h
+	HandleNekMesh1D.h
+	HandleNekMesh2D.h
+	HandleNekMesh3D.h
+	Smoothie.h
+	SmoothieSIAC.h
+	SmoothieSIAC1D.h
+	SmoothieSIAC2D.h
+	SmoothieSIAC3D.h
+	element_centre_size.h
+
+)
+
+SET(LSIACSources
+	NektarBaseClass.cpp
+	GeneralBSplines.cpp
+	CentralBSplines.cpp
+	SIACFilter.cpp
+	SymmetricSIAC.cpp
+	OneSidedSIAC.cpp
+	HandleNekMesh.cpp
+	HandleNekMesh1D.cpp
+	HandleNekMesh2D.cpp
+	HandleNekMesh3D.cpp
+	SmoothieSIAC.cpp
+	SmoothieSIAC1D.cpp
+	SmoothieSIAC2D.cpp
+	SmoothieSIAC3D.cpp
+)
+
+ADD_NEKTAR_LIBRARY(LSIAC lib ${NEKTAR_LIBRARY_TYPE} ${LSIACSources} ${LSIACHeaders})
+
+ADD_DEFINITIONS(-DLSIAC_EXPORTS)
+
+TARGET_LINK_LIBRARIES(LSIAC LINK_PUBLIC Collections)
+
+INSTALL(DIRECTORY ./ DESTINATION ${NEKTAR_INCLUDE_DIR}/LSIAC COMPONENT dev FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

library/LSIAC/CentralBSplines.cpp

+#include "CentralBSplines.h"
+
+namespace Nektar
+{
+    namespace LSIAC
+    {
+        CentralBSplines::CentralBSplines(int Order, NekDouble Shift, NekDouble scale)
+            : GeneralBSplines(Order), m_shift(Shift), m_scale(scale)
+        {
+            Array<OneD,NekDouble> knots(Order+1,0.0);
+            for (int i =0; i < Order+1;i++)
+                {
+                    knots[i] = -Order/2.0 + i;
+                }
+            this->SetKnotVector(knots);
+            GeneralBSplines::SetKnotVector(knots);
+        }
+
+        CentralBSplines::CentralBSplines(int Order)
+            : GeneralBSplines(Order), m_shift(0.0), m_scale(1.0)
+        {
+            Array<OneD,NekDouble> knots(Order+1,0.0);
+            for (int i =0; i < Order+1;i++)
+                {
+                    knots[i] = -Order/2.0 + i;
+                }
+            this->SetKnotVector(knots);
+            GeneralBSplines::SetKnotVector(knots);
+        }
+
+        CentralBSplines::CentralBSplines(int Order, NekDouble shift)
+            : GeneralBSplines(Order), m_shift(shift), m_scale(1.0)
+        {
+            Array<OneD,NekDouble> knots(Order+1,0.0);
+            for (int i =0; i < Order+1;i++)
+                {
+                    knots[i] = -Order/2.0 + i;
+                }
+            this->SetKnotVector(knots);
+            GeneralBSplines::SetKnotVector(knots);
+        }
+
+        bool CentralBSplines::EvaluateBSplines(const Array<OneD,NekDouble> &t_pos,
+                                               Array<OneD,NekDouble> &t_val, const NekDouble shift,
+                                               const NekDouble meshScaling) const
+        {
+            //For central BSplines. The j is always zero.
+            //The spline needs to be shifted and evaluated at x_pos locations.
+            int nq = t_pos.num_elements();
+            GeneralBSplines::EvaluateBSplines(t_pos, 0,t_val,shift, meshScaling);
+            return true;
+        }
+
+    }
+}

library/LSIAC/CentralBSplines.h

+#ifndef LSIAC_CENTRALBSPLINES_H
+#define LSIAC_CENTRALBSPLINES_H
+
+#include "GeneralBSplines.h"
+
+/// The class evaluates CentralB-Splines at given locations.
+/**	This calss automatically calculates knot positions, when using the order specified.
+	All the knot positions calcualted will range -(k+1)/2 to k+1)/2
+*/
+namespace Nektar
+{
+namespace LSIAC
+{
+class CentralBSplines: public GeneralBSplines{
+
+	private:
+	protected:
+	public:
+		NekDouble m_shift;
+		NekDouble m_scale;
+
+	CentralBSplines( int Order );
+	CentralBSplines( int Order, NekDouble shift);
+	CentralBSplines( int Order, NekDouble shift, NekDouble scale);
+
+	bool GetShift( NekDouble &shift) const;
+	bool SetShift( NekDouble shift);
+
+	bool GetScale( NekDouble &scale) const;
+	bool SetScale( NekDouble scale);
+
+	bool SetOrder( int Order);
+	int  GetOrder( ) const;
+	
+//	bool EvaluateBSplines( const Array<OneD,NekDouble> &t_pos,const NekDouble shift, Array<OneD,NekDouble> &t_vals)const;
+	bool EvaluateBSplines( const Array<OneD,NekDouble> &t_pos, Array<OneD,NekDouble> &t_vals, 
+							const NekDouble shift=0.0, const NekDouble meshScaling=1.0 )const;
+//	bool EvaluateBSplines( Array<OneD,NekDouble> &t_pos, Array<OneD,NekDouble> &t_vals);
+
+};
+
+
+}
+}
+
+#endif

library/LSIAC/FilterType.h

+#ifndef FILTERTYPE_H
+#define FILTERTYPE_H
+
+namespace Nektar
+{
+    namespace LSIAC
+    {
+
+        namespace SIACUtilities
+        {
+            //! This enum specifies the kind of filter to be used by the object.
+            /*!  The user can specify any type of smoothing or derivative filter.
+         *  All the properties of memeber depends on this enum. This enum cannot be changed after the object has been initialized.
+         *  A new enum needs to be added to the list to create any new type of filter.
+         */
+            enum FilterType
+            {
+                eNONE,                   /**< enum value 1 */
+                eSYM_2kp1,               /**< enum value 2, Will apply filter only when symmetric filter is possible. */
+                eSYM_4kp1,
+                eSYM,
+                eSYM_2kp1_1SIDED_2kp1,   /**< enum value 3, Simplest form of SIAC and oneSided filter. */
+                eSYM_2kp1_1SIDED_2kp2,   /**< enum value 4, Xli 2k+1 CentralBspl and 1 GeneralBSlp only at borders appropriately*/
+                eSYM_2kp1_1SIDED_XLI,
+                eSYM_2kp1_1SIDED_4kp1,   /**< enum value 5, SRV's 4k+1 CentralBspl at boundaries.*/
+                eSYM_2kp1_1SIDED_SRV,   /**< enum value 5, SRV's 4k+1 CentralBspl at boundaries.*/
+                eSYM_1SIDED_BASIC,
+                eSYM_1SIDED_XLI,
+                eSYM_1SIDED,
+                eSYM_DER_2kp1,
+                eSYM_DER,
+                eSYM_DER_1SIDED,
+                eSYM_DER_2kp1_1SIDED_2kp1,   /**< enum value 6, Simple Derivative filter*/
+                eSYM_DER_2kp1_1SIDED_2kp2,   /**< enum value 7, SRV's Derivative filter */
+                eSYM_DER_2kp1_1SIDED_4kp1,    /**< enum value 8, Xli's Derivative fitler */
+                eSYM_DER_1SIDED_SRV,
+                eSYM_DER_1SIDED_XLi,
+            };
+
+            struct SMOOTHIE_Status {
+                int count;
+                int cancelled;
+                int SOURCE;
+                int TAG;
+                int ERROR;
+            };
+        }
+
+
+    }
+}
+
+#endif

library/LSIAC/GeneralBSplines.cpp

+#include "GeneralBSplines.h"
+
+namespace Nektar
+{
+namespace LSIAC
+{
+
+GeneralBSplines::GeneralBSplines( const int order)
+{
+	this->SetOrder( order);
+}
+
+GeneralBSplines::GeneralBSplines( const Array<OneD,NekDouble> &knots, const int order)
+{
+	this->SetKnotVector( knots);
+	this->SetOrder( order);
+}
+
+bool GeneralBSplines::SetKnotVector( const Array<OneD,NekDouble> &knots)
+{
+	this->m_knotVector=knots;
+	return true;
+}
+
+bool GeneralBSplines::GetKnotVector( Array<OneD,NekDouble> &knots)
+{
+	knots = this->m_knotVector;
+	return true;
+}
+
+
+int GeneralBSplines::GetOrder() const
+{
+	return this->m_order;
+}
+
+bool GeneralBSplines::SetOrder(const int order)
+{
+	this->m_order = order;
+	return true;	
+}
+
+bool GeneralBSplines::EvaluateBSplines (const Array<OneD,NekDouble> &t_pos,const std::vector<NekDouble> &kvec,const int j_th, Array<OneD,NekDouble> &t_values, const NekDouble shift,
+			const NekDouble meshScaling)const
+{
+	return BSplinesBasis( t_pos, kvec, m_order-1, j_th, t_values, shift, meshScaling);
+}	
+
+bool GeneralBSplines::EvaluateBSplines( const Array<OneD,NekDouble> &t_pos, const int j_th,Array<OneD,NekDouble> &t_values,const NekDouble shift, const NekDouble meshScaling) const
+{
+	
+	BSplinesBasis( t_pos, m_order-1, j_th, t_values,shift,meshScaling);
+	return true;	
+}
+
+
+bool GeneralBSplines::EvaluateBSplines( const Array<OneD,NekDouble> &t_pos,const Array<OneD,NekDouble> &knots,int j_th,Array<OneD,NekDouble> &t_values)
+{
+	cout << "This is stub. Need to be coded." << endl;	
+	return true;
+}
+
+bool GeneralBSplines::BSplinesBasis( const Array<OneD,NekDouble> &t_pos, const int k,const int j, 
+					Array<OneD,NekDouble> &t_val, const NekDouble shift,
+					const NekDouble meshScaling)const {
+	// Note here Order of BSplines are k+1. 
+	// This is done to follow the paper. Sorry if causes confusion.
+        if (0==k) {
+		if ( j>=0 && m_knotVector.num_elements()-1){
+			for (int i =0; i < t_pos.num_elements(); i++)
+			{
+				//if( (t_pos[i]/meshScaling-shift >= m_knotVector[j]) & (t_pos[i]/meshScaling-shift < m_knotVector[j+1]))
+				if( ( (t_pos[i]-shift)/meshScaling>= m_knotVector[j]) & ( (t_pos[i]-shift)/meshScaling< m_knotVector[j+1]))
+				{
+					t_val[i] = 	1.0;
+				}else{
+					t_val[i] = 0.0;
+				}
+			}
+		}else{
+			//Add Assert here. Should not come here.
+			printf("Add Asset here");
+		}
+	}else{
+		NekDouble x_eval, w_jlt,w_jlt1;
+		Array<OneD,NekDouble> Bspl_k1_j, Bspl_k1_j1;
+		Bspl_k1_j = Array<OneD,NekDouble>(t_pos.num_elements(),0.0);
+		Bspl_k1_j1 = Array<OneD,NekDouble>(t_pos.num_elements(),0.0);
+		this->BSplinesBasis( t_pos, k-1, j, Bspl_k1_j,shift,meshScaling);
+		this->BSplinesBasis( t_pos, k-1, j+1, Bspl_k1_j1,shift,meshScaling);
+		for (int i =0; i < t_pos.num_elements(); i++)
+		{
+			//x_eval = t_pos[i]/meshScaling-shift;
+			x_eval = (t_pos[i]-shift)/meshScaling;
+                        if( abs(m_knotVector[j+k] - m_knotVector[j]) <= 1e-8){
+				w_jlt =0.0;
+			}else{
+				w_jlt = ( x_eval - m_knotVector[j] )/(m_knotVector[j+k]-m_knotVector[j]);
+			}
+			if( abs(m_knotVector[j+k+1] - m_knotVector[j+1]) <= 1e-8){
+				w_jlt1 =0.0;
+			}else{
+				w_jlt1 = ( x_eval - m_knotVector[j+1] )/(m_knotVector[j+k+1]-m_knotVector[j+1]);
+			}
+			t_val[i] = w_jlt*Bspl_k1_j[i] + (1-w_jlt1)*Bspl_k1_j1[i];
+                }
+	}
+	return true;
+}
+
+bool GeneralBSplines::BSplinesBasis( const Array<OneD,NekDouble> &t_pos,const vector<NekDouble> &kVec,
+					const int k,const int j, Array<OneD,NekDouble> &t_val, const NekDouble shift,
+					const NekDouble meshScaling)const {
+	// Note here Order of BSplines are k+1. 
+	// This is done to follow the paper. Sorry if causes confusion.
+        if (0==k) {
+                if ( j>=0 && kVec.size()-1){
+			for (int i =0; i < t_pos.num_elements(); i++)
+			{
+                                if( ( (t_pos[i]-shift)/meshScaling>= kVec[j]) & ( (t_pos[i]-shift)/meshScaling< kVec[j+1]))
+				{
+					t_val[i] = 	1.0;
+                                }else{
+					t_val[i] = 0.0;
+                                }
+			}
+		}else{
+			//Add Assert here. Should not come here.
+			printf("Add Asset here");
+		}
+	}else{
+		NekDouble x_eval, w_jlt,w_jlt1;
+		Array<OneD,NekDouble> Bspl_k1_j, Bspl_k1_j1;
+		Bspl_k1_j = Array<OneD,NekDouble>(t_pos.num_elements(),0.0);
+		Bspl_k1_j1 = Array<OneD,NekDouble>(t_pos.num_elements(),0.0);
+                this->BSplinesBasis( t_pos,kVec, k-1, j, Bspl_k1_j,shift,meshScaling);
+                this->BSplinesBasis( t_pos, kVec, k-1, j+1, Bspl_k1_j1,shift,meshScaling);
+		for (int i =0; i < t_pos.num_elements(); i++)
+		{
+			//x_eval = t_pos[i]/meshScaling-shift;
+			x_eval = (t_pos[i]-shift)/meshScaling;
+			if( abs(kVec[j+k] - kVec[j]) <= 1e-8){
+				w_jlt =0.0;
+			}else{
+				w_jlt = ( x_eval - kVec[j] )/(kVec[j+k]-kVec[j]);
+			}
+			if( abs(kVec[j+k+1] - kVec[j+1]) <= 1e-8){
+				w_jlt1 =0.0;
+			}else{
+				w_jlt1 = ( x_eval - kVec[j+1] )/(kVec[j+k+1]-kVec[j+1]);
+			}
+			t_val[i] = w_jlt*Bspl_k1_j[i] + (1-w_jlt1)*Bspl_k1_j1[i];
+                }
+	}
+	return true;
+}
+
+
+}
+}

library/LSIAC/GeneralBSplines.h

+#ifndef LSIAC_GENERALBSPLINES_H
+#define LSIAC_GENERALBSPLINES_H
+
+
+#include <iostream>
+#include <vector>
+#include "BSplines.h"
+
+
+/// This class evaluates any general BSplines at given location when knots and order are specified.
+
+/** To evaluate General BSplines at any location, one needs knots and Order of Bsplines.
+*/
+
+namespace Nektar
+{
+namespace LSIAC
+{
+
+class GeneralBSplines: public BSplines{
+    // data
+	public:
+		Array<OneD,NekDouble> m_knotVector;
+		int m_order;
+	// functions
+	private:
+		bool BSplinesBasis(const Array<OneD,NekDouble> &t_pos, const int k,const int j, 
+					Array<OneD,NekDouble> &t_val, const NekDouble scaling=0.0, const NekDouble meshScaling=1.0)const;
+
+		bool BSplinesBasis(const Array<OneD,NekDouble> &t_pos,const vector<NekDouble> &kvector, const int k,const int j, 
+					Array<OneD,NekDouble> &t_val, const NekDouble scaling=0.0, const NekDouble meshScaling=1.0)const;
+	protected:
+	public:
+	GeneralBSplines(const int order);
+	GeneralBSplines(const Array<OneD,NekDouble> &knots,const int order);
+
+	bool SetKnotVector(const Array<OneD,NekDouble> &knots);
+
+	bool GetKnotVector ( Array<OneD,NekDouble> &knots);
+	
+	bool SetOrder(const int order);
+
+	int  GetOrder() const;
+
+	//Not used yet!
+	bool EvaluateBSplines(const  Array<OneD,NekDouble> &t_pos,const Array<OneD,NekDouble> &knots,const int j_th,
+					Array<OneD,NekDouble> &t_values);
+
+	//Not used yet!
+	bool EvaluateBSplines(const  Array<OneD,NekDouble> &t_pos,const Array<OneD,NekDouble> &knots,const int j_th, 
+					const NekDouble shift, Array<OneD,NekDouble> &t_values);
+
+	bool EvaluateBSplines (const Array<OneD,NekDouble> &t_pos,const int j_th, Array<OneD,NekDouble> &t_values,
+								const NekDouble shift=0.0, const NekDouble meshScaling=1.0)const;	
+
+	bool EvaluateBSplines (const Array<OneD,NekDouble> &t_pos,const std::vector<NekDouble> &kvec, const int j_th,
+			Array<OneD,NekDouble> &t_values, const NekDouble shift=0.0, const NekDouble meshScaling=1.0)const;	
+
+};
+
+}
+}
+
+#endif

library/LSIAC/HandleMesh.h

+#ifndef HANDLEMESH_H
+#define HANDLEMESH_H
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include "NektarBaseClass.h"
+using namespace std;
+
+
+/// This is base class to handle different type of meshes.
+
+/** Each subclass can handle one particular format of the mesh.
+        This class allows users to add their own mesh format in the future purpose.
+*/
+
+namespace Nektar
+{
+    namespace LSIAC
+    {
+
+        class HandleMesh: public NektarBaseClass{
+        public:
+        protected:
+
+
+
+
+
+
+
+
+
+
+
+        };
+
+    }
+}
+
+#endif

library/LSIAC/HandleNekMesh.cpp

+#include "HandleNekMesh.h"
+
+namespace Nektar
+	 {
+			namespace LSIAC
+				 {
+						bool HandleNekMesh::EvaluateAt(const Array<OneD,Array<OneD,NekDouble> >point,
+																					 const int gID, const int eID, Array<OneD,NekDouble> &values,int varNum)
+							 {
+
+									if (!m_expansions[varNum]->GetPhysState