Commit 72cad2e3 authored by Julia Docampo Sanchez's avatar Julia Docampo Sanchez

fixed to work with latest nektar version

parent 49ab3b77
# Main library sub-directories, required by all of Nektar++.
SUBDIRS(GlobalMapping LibUtilities LocalRegions Collections MultiRegions
SpatialDomains StdRegions)
SpatialDomains StdRegions LSIAC)
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/library)
IF (NEKTAR_BUILD_UNIT_TESTS)
......
#ifndef LSIAC_BSPLINES_H
#define LSIAC_BSPLINES_H
#include <iostream>
#include <vector>
#include "LSIACPostProcessor.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 LSIACPostProcessor{
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
SET(LSIAC_SOURCES
fieldDerivatives.cpp
LSIACPostProcessor.cpp
GeneralBSplines.cpp
CentralBSplines.cpp
SIACKernel.cpp
SymmetricKernel.cpp
OneSidedKernel.cpp
HandleNekMesh.cpp
HandleNekMesh1D.cpp
HandleNekMesh2D.cpp
HandleNekMesh3D.cpp
SIACFilter.cpp
SIACFilter1D.cpp
SIACFilter2D.cpp
SIACFilter3D.cpp
ppio.cpp
)
SET(LSIAC_HEADERS
LSIACPostProcessor.h
SIACenumerators.h
GeneralBSplines.h
CentralBSplines.h
SIACKernel.h
SymmetricKernel.h
OneSidedKernel.h
HandleMesh.h
HandleNekMesh.h
HandleNekMesh1D.h
HandleNekMesh2D.h
HandleNekMesh3D.h
SIACpostprocessor.h
SIACFilter.h
SIACFilter1D.h
SIACFilter2D.h
SIACFilter3D.h
)
ADD_DEFINITIONS(-DLSIAC_EXPORTS)
ADD_NEKTAR_LIBRARY(LSIAC
SOURCES ${LSIAC_SOURCES}
HEADERS ${LSIAC_HEADERS}
SUMMARY "LSIAC Filtering library"
DESCRIPTION "This library implements the LSIAC postprocessor for 2D and 3D fields. Contains the kernel class aswell as the mesh interaction functions."
)
#TARGET_LINK_LIBRARIES(LSIAC LINK_PUBLIC Collections)
INSTALL(DIRECTORY ./ DESTINATION ${NEKTAR_INCLUDE_DIR}/LSIAC COMPONENT dev FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp")
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE QtCreatorProject>
<!-- Written by QtCreator 3.0.1, 2017-07-05T11:26:18. -->
<qtcreator>
<data>
<variable>ProjectExplorer.Project.ActiveTarget</variable>
<value type="int">-1</value>
</data>
<data>
<variable>ProjectExplorer.Project.EditorSettings</variable>
<valuemap type="QVariantMap">
<value type="bool" key="EditorConfiguration.AutoIndent">true</value>
<value type="bool" key="EditorConfiguration.AutoSpacesForTabs">false</value>
<value type="bool" key="EditorConfiguration.CamelCaseNavigation">true</value>
<valuemap type="QVariantMap" key="EditorConfiguration.CodeStyle.0">
<value type="QString" key="language">Cpp</value>
<valuemap type="QVariantMap" key="value">
<value type="QByteArray" key="CurrentPreferences">CppGlobal</value>
</valuemap>
</valuemap>
<valuemap type="QVariantMap" key="EditorConfiguration.CodeStyle.1">
<value type="QString" key="language">QmlJS</value>
<valuemap type="QVariantMap" key="value">
<value type="QByteArray" key="CurrentPreferences">QmlJSGlobal</value>
</valuemap>
</valuemap>
<value type="int" key="EditorConfiguration.CodeStyle.Count">2</value>
<value type="QByteArray" key="EditorConfiguration.Codec">UTF-8</value>
<value type="bool" key="EditorConfiguration.ConstrainTooltips">false</value>
<value type="int" key="EditorConfiguration.IndentSize">4</value>
<value type="bool" key="EditorConfiguration.KeyboardTooltips">false</value>
<value type="bool" key="EditorConfiguration.MouseNavigation">true</value>
<value type="int" key="EditorConfiguration.PaddingMode">1</value>
<value type="bool" key="EditorConfiguration.ScrollWheelZooming">true</value>
<value type="int" key="EditorConfiguration.SmartBackspaceBehavior">1</value>
<value type="bool" key="EditorConfiguration.SpacesForTabs">true</value>
<value type="int" key="EditorConfiguration.TabKeyBehavior">2</value>
<value type="int" key="EditorConfiguration.TabSize">8</value>
<value type="bool" key="EditorConfiguration.UseGlobal">true</value>
<value type="int" key="EditorConfiguration.Utf8BomBehavior">1</value>
<value type="bool" key="EditorConfiguration.addFinalNewLine">true</value>
<value type="bool" key="EditorConfiguration.cleanIndentation">true</value>
<value type="bool" key="EditorConfiguration.cleanWhitespace">true</value>
<value type="bool" key="EditorConfiguration.inEntireDocument">false</value>
</valuemap>
</data>
<data>
<variable>ProjectExplorer.Project.PluginSettings</variable>
<valuemap type="QVariantMap"/>
</data>
<data>
<variable>ProjectExplorer.Project.TargetCount</variable>
<value type="int">0</value>
</data>
<data>
<variable>ProjectExplorer.Project.Updater.EnvironmentId</variable>
<value type="QByteArray">{a578d475-18b0-40f9-843f-5eee755f73a3}</value>
</data>
<data>
<variable>ProjectExplorer.Project.Updater.FileVersion</variable>
<value type="int">15</value>
</data>
</qtcreator>
#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;
}
}
}
#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
#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;
}
}
}
#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
#ifndef HANDLEMESH_H
#define HANDLEMESH_H
#include <iostream>
#include <vector>
#include <string>
#include "LSIACPostProcessor.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 LSIACPostProcessor{
public:
protected:
};
}
}
#endif
This diff is collapsed.
This diff is collapsed.
#include "HandleNekMesh1D.h"
#include <MultiRegions/DisContField1D.h>
#include <MultiRegions/ContField1D.h>
#include <SpatialDomains/MeshGraph1D.h>
#include <cmath>
/*
HandleNekMesh1D::HandleNekMesh1D (const string meshFileName, const string fldFileName):
HandleNekMesh(meshFileName,fldFileName)
{
//LibUtilitites::SessionReader::CreateInstance
}
*/
//!This function given range of tmin and tmax returns the element break points.
/*
\param xcen_offset,ycen_offset,zcen_offset,direction,t_offsetmin t_offset_max
\param [out] xPos,yPos,zPos,tPos
This function makes few assumptions for simplicity.
-> This function does not gaurentee if all of breakpoints tmin and tmax are returned.
-> This function always includes tmin and tmax as break points while returning.
-> There is no significance to bool in this function. It always return true.
-> This function assumes all the seg element are non-curve elements.
*/
namespace Nektar
{
namespace LSIAC
{
bool HandleNekMesh1D::v_LoadMesh(string var)
{
SpatialDomains::ExpansionMap expansions = m_graph->GetExpansions();
// cout << "expansion size: " <<expansions.size() << endl;
m_expansions.push_back(MemoryManager<MultiRegions::DisContField1D>
::AllocateSharedPtr(m_session, m_graph,var));
return true;
}
bool HandleNekMesh1D::v_LoadData( string filename, vector<string> &variables )
{
SpatialDomains::ExpansionMap expansions = m_graph->GetExpansions();
// cout << "expansion size: " <<expansions.size() << endl;
for (int i =0; i < variables.size(); i++)
{
m_expansions.push_back(MemoryManager<MultiRegions::DisContField1D>
::AllocateSharedPtr(m_session, m_graph,variables[i] ));
}
std::vector<LibUtilities::FieldDefinitionsSharedPtr> rFieldDef;
std::vector<std::vector<NekDouble> > rFieldData;
Array<OneD,int> ElementGIDs(expansions.size());
SpatialDomains::ExpansionMap::const_iterator expIt;
int i=0;
for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
{
ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID();
}
Import( filename, rFieldDef, rFieldData,
LibUtilities::NullFieldMetaDataMap,ElementGIDs);
for(int i =0 ; i < rFieldDef.size() ; i++)
{
for (int e =0; e< variables.size(); e++)
{
m_expansions[e]->ExtractDataToCoeffs(rFieldDef[i], rFieldData[i],
variables[e], m_expansions[e]->UpdateCoeffs());
}
}
for ( int i =0; i < m_expansions.size(); i++)
{
m_expansions[i]->BwdTrans ( m_expansions[i]->GetCoeffs(),
m_expansions[i]->UpdatePhys());
}
return true;
}
// This function does not verify if tmin to tmax is properly given.
bool HandleNekMesh1D::v_GetBreakPts( const NekDouble xcen_offset,const NekDouble ycen_offset,
const NekDouble zcen_offset, Array<OneD,NekDouble> &direction, const NekDouble tmin,