Commit 7097a4ea authored by Ruben Undheim's avatar Ruben Undheim

New upstream version 0.0.35+ds.1

parent aa7abb5c
......@@ -6,7 +6,7 @@ ELSE()
SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" )
ENDIF()
PROJECT( AppCSXCAD CXX)
PROJECT( AppCSXCAD CXX C)
cmake_minimum_required(VERSION 2.8)
......@@ -19,7 +19,7 @@ IF(EXISTS ${PROJECT_SOURCE_DIR}/localConfig.cmake)
ENDIF()
# default
set(VERSION "v0.2.1")
set(VERSION "v0.2.2")
# add git revision
IF(EXISTS ${PROJECT_SOURCE_DIR}/.git )
......@@ -95,6 +95,10 @@ find_path(QCSXCAD_INCLUDE_DIR
message(STATUS "QCSXCAD_INCLUDE_DIR: ${QCSXCAD_INCLUDE_DIR}" )
INCLUDE_DIRECTORIES( ${QCSXCAD_INCLUDE_DIR} )
# hdf5
find_package(HDF5 1.8 COMPONENTS C HL REQUIRED)
INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS})
# vtk
if (WIN32)
find_package(VTK 6.1 REQUIRED)
......@@ -149,6 +153,8 @@ endif()
TARGET_LINK_LIBRARIES( AppCSXCAD
${CSXCAD_LIBRARIES}
${QCSXCAD_LIBRARIES}
${HDF5_LIBRARIES}
${HDF5_HL_LIBRARIES}
${QT_LIBRARIES}
${vtk_LIBS}
)
......
......@@ -6,14 +6,14 @@ ELSE()
SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" )
ENDIF()
PROJECT(CSXCAD CXX)
PROJECT(CSXCAD CXX C)
cmake_minimum_required(VERSION 2.8)
# default
set(LIB_VERSION_MAJOR 0)
set(LIB_VERSION_MINOR 6)
set(LIB_VERSION_PATCH 1)
set(LIB_VERSION_PATCH 2)
set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH})
set(VERSION "v${LIB_VERSION_STRING}")
......@@ -91,7 +91,7 @@ find_package(TinyXML REQUIRED)
ADD_DEFINITIONS( -DTIXML_USE_STL )
find_package(HDF5 1.8 COMPONENTS C HL REQUIRED)
INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIR})
INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS})
link_directories(${HDF5_LIBRARY_DIRS})
# hdf5 compat
ADD_DEFINITIONS( -DH5_USE_16_API )
......
......@@ -89,6 +89,7 @@ TARGET_LINK_LIBRARIES( CSXCAD
${fparser_LIBRARIES}
${TinyXML_LIBRARIES}
${HDF5_LIBRARIES}
${HDF5_HL_LIBRARIES}
CGAL
${Boost_LIBRARIES}
${vtk_LIBS}
......
......@@ -56,9 +56,6 @@ CSPrimBox::~CSPrimBox()
bool CSPrimBox::GetBoundBox(double dBoundBox[6], bool PreserveOrientation)
{
// if ( (m_MeshType!=m_PrimCoordSystem) && (m_PrimCoordSystem!=UNDEFINED_CS))
// std::cerr << "GetBoundBox::GetBoundBox: Warning: The bounding box for this object is not calculated properly... " << std::endl;
const double* start = m_Coords[0].GetCoords(m_MeshType);
const double* stop = m_Coords[1].GetCoords(m_MeshType);
......@@ -81,6 +78,9 @@ bool CSPrimBox::GetBoundBox(double dBoundBox[6], bool PreserveOrientation)
dBoundBox[2*i]=dBoundBox[2*i+1];
dBoundBox[2*i+1]=help;
}
if ( (m_MeshType!=m_PrimCoordSystem) && (m_PrimCoordSystem!=UNDEFINED_CS))
// if the box is defined in a coordinate system other than the expected one, this BB is invalid
return false;
return true;
}
......
......@@ -123,6 +123,13 @@ void CSPrimitives::Init()
m_BoundBoxValid = false;
}
CSTransform* CSPrimitives::GetTransform()
{
if (m_Transform==NULL)
m_Transform = new CSTransform(clParaSet);
return m_Transform;
}
void CSPrimitives::SetProperty(CSProperties *prop)
{
if ((clProperty!=NULL) && (clProperty!=prop))
......@@ -146,8 +153,9 @@ int CSPrimitives::IsInsideBox(const double *boundbox)
return 0; // unable to decide with an invalid bounding box
if ((this->GetBoundBoxCoordSystem()!=UNDEFINED_CS) && (this->GetBoundBoxCoordSystem()!=this->GetCoordInputType()))
return 0; // unable to decide if coordinate system do not match
if (this->GetTransform()!=NULL)
return 0; // unable to decide if a transformation is used
if (m_Transform!=NULL)
if (m_Transform->HasTransform())
return 0; // unable to decide if a transformation is used
for (int i=0;i<3;++i)
{
......
......@@ -177,7 +177,8 @@ public:
//! Read the coordinate system for this primitive (may be different to the input mesh type) \sa GetCoordInputType
CoordinateSystem GetCoordinateSystem() const {return m_PrimCoordSystem;}
CSTransform* GetTransform() const {return m_Transform;}
//! Get the CSTransform if it exists already or create a new one
CSTransform* GetTransform();
//! Show status of this primitve
virtual void ShowPrimitiveStatus(std::ostream& stream);
......
......@@ -73,6 +73,11 @@ void CSTransform::Reset()
MakeUnitMatrix(m_Inv_TMatrix);
}
bool CSTransform::HasTransform()
{
return (m_TransformList.size()>0);
}
void CSTransform::Invert()
{
//make sure the inverse matrix is up to date...
......
......@@ -85,6 +85,9 @@ public:
void Reset();
//! Check if this CSTransform has any transformations
bool HasTransform();
//! All subsequent operations will be occur before the previous operations (not the default).
void SetPreMultiply() {m_PostMultiply=false;}
//! All subsequent operations will be after the previous operations (default).
......
......@@ -22,7 +22,7 @@ ENDIF()
# default
set(LIB_VERSION_MAJOR 0)
set(LIB_VERSION_MINOR 6)
set(LIB_VERSION_PATCH 1)
set(LIB_VERSION_PATCH 2)
set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH})
set(VERSION "v${LIB_VERSION_STRING}")
......@@ -107,7 +107,7 @@ endif()
message(STATUS "Found package VTK. Using version " ${VTK_VERSION})
include(${VTK_USE_FILE})
INCLUDE_DIRECTORIES (${VTK_INCLUDE_DIR})
INCLUDE_DIRECTORIES (${VTK_INCLUDE_DIRS})
# Qt
SET(RESOURCES resources.qrc)
......
......@@ -778,7 +778,7 @@ vtkActor* VTKPrimitives::AddPolyData(vtkPolyData* polydata, double *dRGB, double
filter->SetTransform(vtrans);
#if VTK_MAJOR_VERSION>=6
m_PolyDataCollection->AddInputData(filter->GetOutput());
m_PolyDataCollection->AddInputConnection(filter->GetOutputPort());
#else
m_PolyDataCollection->AddInput(filter->GetOutput());
#endif
......@@ -812,7 +812,7 @@ vtkActor* VTKPrimitives::AddPolyData(vtkAlgorithmOutput* polydata_port, double *
filter->SetTransform(vtrans);
#if VTK_MAJOR_VERSION>=6
m_PolyDataCollection->AddInputData(filter->GetOutput());
m_PolyDataCollection->AddInputConnection(filter->GetOutputPort());
#else
m_PolyDataCollection->AddInput(filter->GetOutput());
#endif
......@@ -855,7 +855,7 @@ void VTKPrimitives::WritePolyData2File(const char* filename, double scale)
if (scale==1.0)
{
#if VTK_MAJOR_VERSION>=6
writer->SetInputData(m_PolyDataCollection->GetOutput());
writer->SetInputConnection(m_PolyDataCollection->GetOutputPort());
#else
writer->SetInput(m_PolyDataCollection->GetOutput());
#endif
......@@ -867,7 +867,7 @@ void VTKPrimitives::WritePolyData2File(const char* filename, double scale)
vtkTransformPolyDataFilter *transformFilter = vtkTransformPolyDataFilter::New();
#if VTK_MAJOR_VERSION>=6
transformFilter->SetInputData(m_PolyDataCollection->GetOutput());
transformFilter->SetInputConnection(m_PolyDataCollection->GetOutputPort());
#else
transformFilter->SetInput(m_PolyDataCollection->GetOutput());
#endif
......@@ -875,7 +875,7 @@ void VTKPrimitives::WritePolyData2File(const char* filename, double scale)
transformFilter->SetTransform(transform);
#if VTK_MAJOR_VERSION>=6
writer->SetInputData(transformFilter->GetOutput());
writer->SetInputConnection(transformFilter->GetOutputPort());
#else
writer->SetInput(transformFilter->GetOutput());
#endif
......@@ -895,7 +895,7 @@ void VTKPrimitives::WritePolyData2STL(const char* filename, double scale)
vtkTriangleFilter* filter = vtkTriangleFilter::New();
#if VTK_MAJOR_VERSION>=6
filter->SetInputData(m_PolyDataCollection->GetOutput());
filter->SetInputConnection(m_PolyDataCollection->GetOutputPort());
#else
filter->SetInput(m_PolyDataCollection->GetOutput());
#endif
......@@ -917,7 +917,7 @@ void VTKPrimitives::WritePolyData2STL(const char* filename, double scale)
transformFilter->SetTransform(transform);
#if VTK_MAJOR_VERSION>=6
writer->SetInputData(transformFilter->GetOutput());
writer->SetInputConnection(transformFilter->GetOutputPort());
#else
writer->SetInput(transformFilter->GetOutput());
#endif
......@@ -937,7 +937,7 @@ void VTKPrimitives::WritePolyData2PLY(const char* filename, double scale)
vtkTriangleFilter* filter = vtkTriangleFilter::New();
#if VTK_MAJOR_VERSION>=6
filter->SetInputData(m_PolyDataCollection->GetOutput());
filter->SetInputConnection(m_PolyDataCollection->GetOutputPort());
#else
filter->SetInput(m_PolyDataCollection->GetOutput());
#endif
......@@ -960,7 +960,7 @@ void VTKPrimitives::WritePolyData2PLY(const char* filename, double scale)
transformFilter->SetTransform(transform);
#if VTK_MAJOR_VERSION>=6
writer->SetInputData(transformFilter->GetOutput());
writer->SetInputConnection(transformFilter->GetOutputPort());
#else
writer->SetInput(transformFilter->GetOutput());
#endif
......
......@@ -6,13 +6,13 @@ ELSE()
SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" )
ENDIF()
PROJECT(openEMS CXX)
PROJECT(openEMS CXX C)
cmake_minimum_required(VERSION 2.8)
# default
set(LIB_VERSION_MAJOR 0)
set(LIB_VERSION_MINOR 0)
set(LIB_VERSION_PATCH 34)
set(LIB_VERSION_PATCH 35)
set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH})
set(VERSION "v${LIB_VERSION_STRING}")
......@@ -21,7 +21,7 @@ IF(EXISTS ${PROJECT_SOURCE_DIR}/localConfig.cmake)
include(${PROJECT_SOURCE_DIR}/localConfig.cmake)
ENDIF()
set(VERSION "v0.0.34")
set(VERSION "v0.0.35")
# add git revision
IF(EXISTS ${PROJECT_SOURCE_DIR}/.git )
......@@ -109,11 +109,8 @@ ADD_DEFINITIONS( -DTIXML_USE_STL )
# hdf5
find_package(HDF5 1.8 COMPONENTS C HL REQUIRED)
INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIR})
link_directories(${HDF5_LIBRARY_DIRS})
# hdf5 compat
ADD_DEFINITIONS( -DH5_USE_16_API )
INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS})
link_directories(${HDF5_LIBRARIES})
# boost
find_package(Boost 1.46 COMPONENTS
......@@ -158,6 +155,8 @@ set(SOURCES
openems.cpp
)
set(PUB_HEADERS openems.h)
# libs
ADD_SUBDIRECTORY( tools )
ADD_SUBDIRECTORY( FDTD )
......@@ -181,6 +180,7 @@ TARGET_LINK_LIBRARIES( openEMS
${fparser_LIBRARIES}
tinyxml
${HDF5_LIBRARIES}
${HDF5_HL_LIBRARIES}
${Boost_LIBRARIES}
${vtk_LIBS}
${MPI_LIBRARIES}
......@@ -204,5 +204,6 @@ if (UNIX)
PERMISSIONS OWNER_WRITE OWNER_READ GROUP_READ WORLD_READ OWNER_EXECUTE GROUP_EXECUTE WORLD_EXECUTE)
endif()
endif ()
INSTALL(FILES ${PUB_HEADERS} DESTINATION include/openEMS)
INSTALL( DIRECTORY matlab DESTINATION share/openEMS )
# TODO mpi, tarball, debug, release
......@@ -135,22 +135,23 @@ void ProcessModeMatch::InitProcess()
m_ModeDist[n] = Create2DArray<double>(m_numLines);
}
bool dualMesh = m_ModeFieldType==1;
unsigned int pos[3] = {0,0,0};
double discLine[3] = {0,0,0};
double gridDelta = 1; // 1 -> mode-matching function is definied in drawing units...
double var[7];
pos[m_ny] = start[m_ny];
discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],m_dualMesh);
discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],dualMesh);
double norm = 0;
double area = 0;
for (unsigned int posP = 0; posP<m_numLines[0]; ++posP)
{
pos[nP] = start[nP] + posP;
discLine[nP] = Op->GetDiscLine(nP,pos[nP],m_dualMesh);
discLine[nP] = Op->GetDiscLine(nP,pos[nP],dualMesh);
for (unsigned int posPP = 0; posPP<m_numLines[1]; ++posPP)
{
pos[nPP] = start[nPP] + posPP;
discLine[nPP] = Op->GetDiscLine(nPP,pos[nPP],m_dualMesh);
discLine[nPP] = Op->GetDiscLine(nPP,pos[nPP],dualMesh);
var[0] = discLine[0] * gridDelta; // x
var[1] = discLine[1] * gridDelta; // y
......@@ -169,7 +170,7 @@ void ProcessModeMatch::InitProcess()
var[5] = sqrt(pow(discLine[0],2)+pow(discLine[2],2)) * gridDelta; // r
var[6] = asin(1)-atan(var[2]/var[3]); //theta (t)
}
area = Op->GetNodeArea(m_ny,pos,m_dualMesh);
area = Op->GetNodeArea(m_ny,pos,dualMesh);
for (int n=0; n<2; ++n)
{
m_ModeDist[n][posP][posPP] = m_ModeParser[n]->Eval(var); //calc mode template
......@@ -227,6 +228,7 @@ double* ProcessModeMatch::CalcMultipleIntegrals()
double field = 0;
double purity = 0;
double area = 0;
bool dualMesh = m_ModeFieldType==1;
int nP = (m_ny+1)%3;
int nPP = (m_ny+2)%3;
......@@ -242,7 +244,7 @@ double* ProcessModeMatch::CalcMultipleIntegrals()
for (unsigned int posPP = 0; posPP<m_numLines[1]; ++posPP)
{
pos[nPP] = start[nPP] + posPP;
area = Op->GetNodeArea(m_ny,pos,m_dualMesh);
area = Op->GetNodeArea(m_ny,pos,dualMesh);
if (m_ModeFieldType==0)
m_Eng_Interface->GetEField(pos,out);
if (m_ModeFieldType==1)
......
......@@ -137,7 +137,7 @@ unsigned int Excitation::GetMaxExcitationTimestep() const
{
FDTD_FLOAT maxAmp=0;
unsigned int maxStep=0;
for (unsigned int n=1; n<Length+1; ++n)
for (unsigned int n=0; n<Length; ++n)
{
if (fabs(Signal_volt[n])>maxAmp)
{
......@@ -152,7 +152,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
{
if (dT==0) return;
Length = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
Length = (unsigned int)ceil(2.0 * 9.0/(2.0*PI*fc) / dT);
if (Length>(unsigned int)nTS)
{
cerr << "Operator::CalcGaussianPulsExcitation: Requested excitation pusle would be " << Length << " timesteps or " << Length * dT << " s long. Cutting to max number of timesteps!" << endl;
......@@ -160,13 +160,11 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
}
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0;
Signal_curr[0]=0.0;
for (unsigned int n=1; n<Length+1; ++n)
Signal_volt = new FDTD_FLOAT[Length];
Signal_curr = new FDTD_FLOAT[Length];
for (unsigned int n=0; n<Length; ++n)
{
double t = (n-1)*dT;
double t = n*dT;
Signal_volt[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
t += 0.5*dT;
Signal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
......@@ -182,12 +180,12 @@ void Excitation::CalcDiracPulsExcitation()
{
if (dT==0) return;
Length = 1;
Length = 2;
// cerr << "Operator::CalcDiracPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt = new FDTD_FLOAT[Length];
Signal_curr = new FDTD_FLOAT[Length];
Signal_volt[0]=0.0;
Signal_volt[1]=1.0;
Signal_curr[0]=0.0;
......@@ -203,11 +201,11 @@ void Excitation::CalcStepExcitation()
{
if (dT==0) return;
Length = 1;
Length = 2;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt = new FDTD_FLOAT[Length];
Signal_curr = new FDTD_FLOAT[Length];
Signal_volt[0]=1.0;
Signal_volt[1]=1.0;
Signal_curr[0]=1.0;
......@@ -228,10 +226,8 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
// cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt[0]=0.0;
Signal_curr[0]=0.0;
Signal_volt = new FDTD_FLOAT[Length];
Signal_curr = new FDTD_FLOAT[Length];
FunctionParser fParse;
fParse.AddConstant("pi", 3.14159265358979323846);
fParse.AddConstant("e", 2.71828182845904523536);
......@@ -242,9 +238,9 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
exit(1);
}
double vars[1];
for (unsigned int n=1; n<Length+1; ++n)
for (unsigned int n=0; n<Length; ++n)
{
vars[0] = (n-1)*dT;
vars[0] = n*dT;
Signal_volt[n] = fParse.Eval(vars);
vars[0] += 0.5*dT;
Signal_curr[n] = fParse.Eval(vars);
......@@ -260,17 +256,16 @@ void Excitation::CalcSinusExcitation(double f0, int nTS)
if (dT==0) return;
if (nTS<=0) return;
Length = (unsigned int)(2.0/f0/dT);
//cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << Length << " timesteps " << Length*dT << "s" << endl;
Length = (unsigned int)round(2.0/f0/dT);
delete[] Signal_volt;
delete[] Signal_curr;
Signal_volt = new FDTD_FLOAT[Length+1];
Signal_curr = new FDTD_FLOAT[Length+1];
Signal_volt = new FDTD_FLOAT[Length];
Signal_curr = new FDTD_FLOAT[Length];
Signal_volt[0]=0.0;
Signal_curr[0]=0.0;
for (unsigned int n=1; n<Length+1; ++n)
for (unsigned int n=1; n<Length; ++n)
{
double t = (n-1)*dT;
double t = n*dT;
Signal_volt[n] = sin(2.0*PI*f0*t);
t += 0.5*dT;
Signal_curr[n] = sin(2.0*PI*f0*t);
......@@ -286,8 +281,8 @@ void Excitation::DumpVoltageExcite(string filename)
file.open( filename.c_str() );
if (file.fail())
return;
for (unsigned int n=1; n<Length+1; ++n)
file << (n-1)*dT << "\t" << Signal_volt[n] << "\n";
for (unsigned int n=0; n<Length; ++n)
file << n*dT << "\t" << Signal_volt[n] << "\n";
file.close();
}
......@@ -297,8 +292,8 @@ void Excitation::DumpCurrentExcite(string filename)
file.open( filename.c_str() );
if (file.fail())
return;
for (unsigned int n=1; n<Length+1; ++n)
file << (n-1)*dT + 0.5*dT << "\t" << Signal_curr[n] << "\n";
for (unsigned int n=0; n<Length; ++n)
file << n*dT + 0.5*dT << "\t" << Signal_curr[n] << "\n";
file.close();
}
......@@ -54,7 +54,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
exc_pos *= (exc_pos<=(int)length);
exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Volt_dir[n];
pos[0]=m_Op_Exc->Volt_index[0][n];
pos[1]=m_Op_Exc->Volt_index[1][n];
......@@ -71,7 +71,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
exc_pos *= (exc_pos<=(int)length);
exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Volt_dir[n];
pos[0]=m_Op_Exc->Volt_index[0][n];
pos[1]=m_Op_Exc->Volt_index[1][n];
......@@ -87,7 +87,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
exc_pos *= (exc_pos<=(int)length);
exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Volt_dir[n];
pos[0]=m_Op_Exc->Volt_index[0][n];
pos[1]=m_Op_Exc->Volt_index[1][n];
......@@ -124,7 +124,7 @@ void Engine_Ext_Excitation::Apply2Current()
exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
exc_pos *= (exc_pos<=(int)length);
exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Curr_dir[n];
pos[0]=m_Op_Exc->Curr_index[0][n];
pos[1]=m_Op_Exc->Curr_index[1][n];
......@@ -141,7 +141,7 @@ void Engine_Ext_Excitation::Apply2Current()
exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
exc_pos *= (exc_pos<=(int)length);
exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Curr_dir[n];
pos[0]=m_Op_Exc->Curr_index[0][n];
pos[1]=m_Op_Exc->Curr_index[1][n];
......@@ -157,7 +157,7 @@ void Engine_Ext_Excitation::Apply2Current()
exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
exc_pos *= (exc_pos>0);
exc_pos %= p;
exc_pos *= (exc_pos<=(int)length);
exc_pos *= (exc_pos<(int)length);
ny = m_Op_Exc->Curr_dir[n];
pos[0]=m_Op_Exc->Curr_index[0][n];
pos[1]=m_Op_Exc->Curr_index[1][n];
......
......@@ -44,7 +44,7 @@ void Engine_Ext_TFSF::DoPostVoltageUpdates()
{
if ( numTS < n )
m_DelayLookup[n]=0;
else if ((numTS-n > length) && (p==0))
else if ((numTS-n >= length) && (p==0))
m_DelayLookup[n]=0;
else
m_DelayLookup[n] = numTS - n;
......@@ -134,7 +134,7 @@ void Engine_Ext_TFSF::DoPostCurrentUpdates()
{
if ( numTS < n )
m_DelayLookup[n]=0;
else if ((numTS-n > length) && (p==0))
else if ((numTS-n >= length) && (p==0))
m_DelayLookup[n]=0;
else
m_DelayLookup[n] = numTS - n;
......
......@@ -210,7 +210,7 @@ bool openEMS_FDTD_MPI::SetupMPI()
if (numProcs!=m_NumProc)
{
if (m_MyID==0)
cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but only " << m_NumProc << " were found! Exit! " << endl;
cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but " << m_NumProc << " were found! Exit! " << endl;
exit(10);
}
......@@ -256,8 +256,8 @@ bool openEMS_FDTD_MPI::SetupMPI()
grid->AddDiscLine(2, m_Original_Grid->GetLine(2,n) );
m_MPI_Op->SetSplitPos(0,m_SplitNumber[0].at(i));
m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(i));
m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(i));
m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(j));
m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(k));
if (i>0)
m_MPI_Op->SetNeighborDown(0,procTable[i-1][j][k]);
......@@ -398,6 +398,7 @@ bool openEMS_FDTD_MPI::SetupProcessing()
{
//type is integral processing --> disable! Needs to be fixed!
cerr << "openEMS_FDTD_MPI::SetupProcessing(): Warning: Processing: " << proc->GetName() << " occures multiple times and is being deactivated..." << endl;
cerr << "openEMS_FDTD_MPI::SetupProcessing(): Note: Processing: Make sure that there are no splits inside probes or sources." << endl;
deactivate = true;
rename = false;
}
......
......@@ -68,8 +68,12 @@ int main(int argc, char *argv[])
}
int EC = FDTD.ParseFDTDSetup(argv[1]);
if(!EC) {
cerr << "openEMS - ParseFDTDSetup failed." << endl;
exit(1);
}
EC = FDTD.SetupFDTD();
if (EC) return EC;
if (EC) exit(EC);
FDTD.RunFDTD();
#ifdef MPI_SUPPORT
......@@ -77,5 +81,5 @@ int main(int argc, char *argv[])
MPI::Finalize();
#endif
return 0;
exit(0);
}
function [delay, fidelity, nf2ff_out] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_0, f_c, varargin)
% [delay, fidelity] = DelayFidelity(nf2ff, port, path, theta, phi, f_lo, f_hi, varargin)
%
%
% This function calculates the time delay from the source port to the phase center of the antenna and the fidelity.
% The fidelity is the similarity between the excitation pulse and the radiated pulse (normalized scalar product).
% The resolution of the delay will be equal to or better than ((f_0 + f_c)*Oversampling)^-1 when using Gaussian excitation.
% Oversampling is an input parameter to InitFDTD. The rows of delay and fidelity correspond to theta and the columns to phi.
%
% input:
% nf2ff: return value of CreateNF2FFBox.
% port: return value of AddLumpedPort
% path: path of the simulation results.
% weight_theta: weight if the E_theta component
% weight_phi: eight of the E_phi component
% -> with both (possibly complex) parameters any polarization can be examined
% theta: theta values to be simulated
% phi: phi values to be simulated
% f_0: center frequency of SetGaussExcite
% f_c: cutoff frequency of SetGaussExcite
%
% variable input:
% 'Center': phase center of the antenna for CalcNF2FF
% 'Radius': radius for CalcNF2FF
% 'Mode': mode CalcNF2FF
%
% example:
% theta = [-180:10:180] * pi / 180;
% phi = [0, 90] * pi / 180;
% [delay, fidelity] = DelayFidelity2(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1);
% figure
% polar(theta.', delay(:,1) * 3e11); % delay in mm
% figure
% polar(theta', (fidelity(:,1)-0.95)/0.05); % last 5 percent of fidelity
%
% Author: Georg Michel
C0 = 299792458;
center = [0, 0, 0];
radius = 1;
nf2ff_mode = 0;
for n=1:2:numel(varargin)
if (strcmp(varargin{n},'Center')==1);
center = varargin{n+1};
elseif (strcmp(varargin{n},'Radius')==1);
radius = varargin{n+1};
elseif (strcmp(varargin{n},'Mode')==1);
nf2ff_mode = varargin{n+1};
end
end
port_ut = load(fullfile(path, port.U_filename));
port_it = load(fullfile(path, port.I_filename));
dt = port_ut(2,1) - port_ut(1,1);
fftsize = 2^(nextpow2(size(port_ut)(1)) + 1);
df = 1 / (dt * fftsize);
uport = fft(port_ut(:, 2), fftsize)(1:fftsize/2+1);
iport = fft(port_it(:, 2), fftsize)(1:fftsize/2+1);
fport = df * (0:fftsize/2);
f_ind = find(fport > (f_0 - f_c ) & fport < (f_0 + f_c));
disp(["frequencies: ", num2str(numel(f_ind))]);
exc_f = uport.' + iport.' * port.Feed_R; %excitation in freq domain
exc_f(!f_ind) = 0;
exc_f /= sqrt(exc_f * exc_f'); % normalization (transposing also conjugates)
nf2ff = CalcNF2FF(nf2ff, path, fport(f_ind), theta, phi, ...
'Center', center, 'Radius', radius, 'Mode', nf2ff_mode);
radfield = weight_theta * cell2mat(nf2ff.E_theta) + weight_phi * cell2mat(nf2ff.E_phi); % rows: theta(f1), columns: phi(f1), phi(f2), ...phi(fn)
radfield = reshape(radfield, [length(nf2ff.theta), length(nf2ff.phi), length(nf2ff.freq)]);
correction = reshape(exp(-2i*pi*nf2ff.r/C0*nf2ff.freq), 1,1,numel(nf2ff.freq)); %dimensions: theta, phi, frequencies
radfield = radfield./correction; % correct for radius delay
% normalize radfield
radnorm = sqrt(dot(radfield, radfield, 3));
radfield ./= radnorm;
%initialize radiated field in fully populated frequency domain
rad_f = zeros([numel(nf2ff.theta), numel(nf2ff.phi), numel(fport)]);
rad_f(:, :, f_ind) = radfield; % assign selected frequencies
exc_f = reshape(exc_f, [1,1,numel(exc_f)]); %make exc_f confomant with rad_f
cr_f = rad_f .* conj(exc_f); % calculate cross correlation
% calculate the cross correlation in time domain (analytic signal)
cr = ifft(cr_f(:, :, 1:end-1), [], 3) * (numel(fport) -1); % twice the FFT normalization (sqrt^2) because product of two normalized functions
%search for the maxiumum of the envelope
[fidelity, delay_ind] = max(abs(cr), [], 3);
delay = (delay_ind - 1) * dt * 2; % double time step because of single-sided FFT
nf2ff_out = nf2ff; %possibly needed for plotting the far field and other things
disp(["DelayFidelity: delay resolution = ", num2str(dt*2e9), "ns"]);
return;
......@@ -21,7 +21,7 @@ function RunOpenEMS(Sim_Path, Sim_File, opts, Settings)
% --engine=fastest fastest available engine (default)