Commit fb8eead4 authored by Julia Docampo Sanchez's avatar Julia Docampo Sanchez

FILTERING :)

parent 27a6ee0c
......@@ -48,18 +48,16 @@ namespace Nektar
return mcoeffs;
}
bool EvaluateKernel(const Array<OneD,NekDouble> &x_pos,
Array<OneD,NekDouble> &vals, const NekDouble meshScaling=1.0,
const NekDouble meshShift=0.0, const bool evalCoeff=false )
Array<OneD,NekDouble> &vals, const NekDouble meshShift=0.0, const bool evalCoeff=false )
{
return v_EvaluateKernel(x_pos, vals, meshScaling,meshShift,evalCoeff );
return v_EvaluateKernel(x_pos, vals,meshShift,evalCoeff );
}
bool GetFilterRange( NekDouble scaling, NekDouble &tmin, NekDouble &tmax,const NekDouble shift=0.0);
bool GetBreakPts( const NekDouble scaling,vector<NekDouble> &valT,const NekDouble shift=0.);
protected:
virtual bool v_EvaluateKernel(const Array<OneD,NekDouble> &x_pos,
Array<OneD,NekDouble> &vals, const NekDouble meshScaling=1.0,
const NekDouble meshShift=0.0, const bool evalCoeff=false )
Array<OneD,NekDouble> &vals, const NekDouble meshShift=0.0, const bool evalCoeff=false )
{
cout << "Assert. Evaluate kernel should be overwritten by a child class " << endl;
return true;
......
......@@ -6,6 +6,7 @@ namespace Nektar
void check_mesh_intersections(conv_pointsSharedPtr &cp,point_info p0,point_info p1,Array<OneD, MultiRegions::DisContField2DSharedPtr> Exp, int dimMesh)
{
cout<< "IN CHEDCK"<<endl;
Array<OneD,NekDouble> dir(dimMesh);
for(int i=0;i<dimMesh;i++)dir[i]=p1.coords[i]-p0.coords[i];
normalise_vector(dir);
......@@ -182,26 +183,11 @@ namespace Nektar
{
case(eVertex):
{
if(cp->kernel_type == e2D)
{
if(cp->new_vertex(p.v_id))
{
for(int i=0;i<cp->dimMesh;i++) cp->meshVertices[i].push_back(p.coords[i]);
cp->vertex_map.push_back(p.v_id);
}
}
else if(cp->kernel_type == e1D)
{
if(cp->new_point(p.coords))
{
for(int i=0;i<cp->dimMesh;i++) cp->meshPoints[i].push_back(p.coords[i]);
}
}
else
{
cout<<"CP has wrong kernel type. Rename"<<endl;
exit(1);
}
}
break;
case(eEdge):
......
......@@ -123,8 +123,9 @@ namespace Nektar
theta=Array<OneD,NekDouble>(dimMesh,0.0);
theta[0] = atan(delta[1]/delta[0]);
theta[1] = theta[0]+M_PI*0.5;
cout<<" element size "<<delta[0]<<"\t"<<delta[1]<<endl;
scaling[0]= delta[0]/fabs(cos(theta[0]));
scaling[1]= delta[0]/fabs(cos(theta[0]));
scaling[1]= delta[1]/fabs(cos(theta[0]));
return TPIntegral(eval_point,pId,filtered_field);
}
......@@ -180,6 +181,8 @@ namespace Nektar
filtered_field[n]=0.;
// SET SCALING
rs=MemoryManager<reference_system>::AllocateSharedPtr(x_bar,theta);
cout<<" DIRECTION "<<theta[0]<<"\t"<<theta[1]<<endl;
cout<<" scaling "<<scaling[0]<<"\t"<<scaling[1]<<endl;
// check that filter support fits
cout<<" CHECK SUPPORT "<<endl;
Array<OneD,NekDouble> vert2(dimMesh);
......@@ -188,6 +191,7 @@ namespace Nektar
{
for(int i=0;i<dimMesh;i++) vert[i]=support_bounds[v][i]*scaling[i];
rs->change_reference_system(vert,vert2,eR2toR1); // change to DG system and get intersections in R1
cout<<" VERT POINT "<<vert2[0]<<"\t"<<vert2[1]<<endl;
if(Exp[0]->GetExpIndex(vert2,NekConstants::kGeomFactorsTol) == -1)
{
for(int f=0;f<dimField;++f)
......@@ -219,12 +223,13 @@ namespace Nektar
scaled_box_matrix[box][v][i]=box_matrix[box][v][i]*scaling[i];
}
rs->change_reference_system(scaled_box_matrix[box][v],rotated_box_matrix[box][v],eR2toR1);
//box_file<<rotated_box_matrix[box][v][0]<<"\t"<<rotated_box_matrix[box][v][1]<<endl;
cout<<rotated_box_matrix[box][v][0]<<"\t"<<rotated_box_matrix[box][v][1]<<endl;
}
}
//box_file.close();
for(int p=0;p<total_kernel_boxes;p++)
{
cout<<" P "<<p<<endl;
int min=Exp[0]->GetExpIndex(rotated_box_matrix[p][0],NekConstants::kGeomFactorsTol);
int max=min;
for(int i=1;i<verts_per_box;i++)
......@@ -235,8 +240,11 @@ namespace Nektar
}
if(min>0) --min;
if(max<Exp[0]->GetNumElmts()) ++max;
cout<<" CP CONSTRU"<<endl;
conv_pointsSharedPtr cp=MemoryManager<convolution_points>::AllocateSharedPtr(rotated_box_matrix[p],max-min+1,min); // structure containing all intersection points stored as vertex, edge or
cout<<" CP CONSTRU"<<endl;
check_element_vertex (min,max,cp,Exp);
cout<<" OK VERTEX"<<endl;
point_info pd(dimMesh);
point_info pd_next(dimMesh);
for(int i=0;i<dimMesh;i++)
......@@ -244,7 +252,9 @@ namespace Nektar
pd.coords[i]=rotated_box_matrix[p][0][i];
pd.aux_coords[i]=pd.coords[i];
}
cout<<" POINT INFO"<<endl;
get_point_info(pd,pd.coords,Exp);
cout<<" COLLECT"<<endl;
collect_point(cp,pd);
for(int b=1;b<=verts_per_box;b++)
{
......@@ -259,7 +269,7 @@ namespace Nektar
for(int i=0;i<dimMesh;i++) pd.coords[i]=rotated_box_matrix[p][b][i];
get_point_info(pd,pd.coords,Exp);
}
//print_intersection_map (cp,eR1);
print_intersection_map (cp,eR1);
Array<OneD,NekDouble> integral_piece;
integrate_region(cp,integral_piece,scaling);
for(int f=0;f<dimField;++f) filtered_field[f]+=integral_piece[f];
......
......@@ -135,11 +135,13 @@ int main(int argc, char *argv[])
for(int d=0;d<dimMesh;++d)
ppPointsList[i][d] = xc[d][i];
}
double tmin=0;
double tmax=TotMeshPoints;
double tmin=500;
double tmax=tmin+1;//TotMeshPoints;
Array<OneD,Array<OneD,NekDouble> > filtered_sol(TotMeshPoints);
for(int pt= tmin;pt<tmax;++pt)
{
ppPointsList[pt][0]= 0.5;
ppPointsList[pt][1]= 0.5;
filtered_sol[pt]=Array<OneD,NekDouble> (totFields);
bool filter=tpf->FilterPoint(ppPointsList[pt],filtered_sol[pt]);
cout<<boolalpha<<filter<<endl;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment