Commit 3608d367 authored by Michael Turner's avatar Michael Turner

Merge branch 'master' into 'ticket/70-Fix-Hdf5-Homogeneous'

# Conflicts:
#   CHANGELOG.md
parents be98663c 0d1375c3
......@@ -23,6 +23,7 @@ v5.0.0
- Refactor ParseUtils to be more consistent (!843)
- Added support for using the distance to a specific region (e.g. outlet) in the
function definitions for the Absorption Forcing (!769)
- Improve performance of DisContField2D::v_ExtractTracePhys (!824)
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
......@@ -34,7 +35,9 @@ v5.0.0
- Fix issue with reading CCM files due to definition of default arrays
rather than a vector (!797)
- Fix inverted triangles and small memory issue in surface meshing (!798)
- Update for the CAD system, more advance self-healing and analysis (!822)
- Additional curve types in GEO reader: BSpline, Circle, Ellipse (!800)
- Fix default command line argument value (!823)
**FieldConvert**:
- Add input module for Semtex field files (!777)
......@@ -45,9 +48,23 @@ v5.0.0
- Enable output to multiple files (!844)
- Allow using xml file without expansion tag in FieldConvert (!849)
**CompressibleFlowSolver**
- Add 3D regression tests (!567)
**Documentation**:
- Added the developer-guide repository as a submodule (!751)
v4.4.2
------
**Library**
- Fix evaluation of points (e.g. HistoryPoints, Interpolation to pts) close to
the interface of two elements (!836)
- Fix deadlock in Hdf5 with homogeneous expansions (!858)
**NekMesh**
- Fix missing periodic boundary meshing and boundary layer mesh adjustment
configurations in 2D (!859)
v4.4.1
------
**Library**
......@@ -65,7 +82,7 @@ v4.4.1
- Fix deadlock with HDF5 input (!786)
- Fix missing entriess in LibUtilities::kPointsTypeStr (!792)
- Fix compiler warnings with CommDataType (!793)
- Fix ability to set default implementation in Collections and added an option
- Fix ability to set default implementation in Collections and added an option
to set eNoCollections in FieldConvert as default (!789)
- Fix performance issue in ProcessIsoContour in relation to memory consumption
(!821)
......@@ -73,7 +90,11 @@ v4.4.1
- Fix available classes being listed multiple times (!817)
- Fix Intel compiler warnings (!837)
- Fix overwriting and backup of chk/fld files on slow file systes (!741)
- Fix deadlock in Hdf5 with homogeneous expansions (!858)
- Fix DriverAdaptive with second order IMEX (!850)
- Fixed typo in eIMEXGear part (!854)
- Added regression tests for IMEXOrder1, IMEXOrder2, IMEXOrder3, MCNAB,
IMEXGear, CNAB, 2nd order IMEX-DIRK, 3rd order IMEX-DIRK (!854)
- Fix bug due to subtractive cancellation in polylib routines (!778)
**FieldConvert:**
- Fix issue with field ordering in the interppointdatatofld module (!754)
......
......@@ -24,19 +24,28 @@ if(NOT DEFINED OCE_DIR)
# Check for OSX needs to come first because UNIX evaluates to true on OSX
if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
if(DEFINED MACPORTS_PREFIX)
find_package(OCE 0.15 QUIET HINTS ${MACPORTS_PREFIX}/Library/Frameworks)
find_package(OCE 0.17 QUIET HINTS ${MACPORTS_PREFIX}/Library/Frameworks)
elseif(DEFINED HOMEBREW_PREFIX)
find_package(OCE 0.15 QUIET HINTS ${HOMEBREW_PREFIX}/Cellar/oce/*)
find_package(OCE 0.17 QUIET HINTS ${HOMEBREW_PREFIX}/Cellar/oce/*)
endif()
elseif(UNIX)
set(OCE_DIR "/usr/local/share/cmake/")
endif()
endif()
find_package(OCE 0.15 QUIET)
find_package(OCE 0.17 QUIET)
if(OCE_FOUND)
message(STATUS "-- OpenCASCADE Community Edition has been found.")
set(OCC_INCLUDE_DIR ${OCE_INCLUDE_DIRS})
message(STATUS "OpenCASCADE Community Edition has been found.")
#check that the OCE package has CAF modules
FIND_LIBRARY(OCC_CAF_LIBRARY TKXCAF ${OCE_INCLUDE_DIRS}/../../lib )
if(OCC_CAF_LIBRARY)
set(OCC_INCLUDE_DIR ${OCE_INCLUDE_DIRS})
else()
message(STATUS "-- OCE does not have CAF libraries, will build from source.")
endif()
else(OCE_FOUND) #look for OpenCASCADE
FIND_PATH(OCC_INCLUDE_DIR Standard_Version.hxx
......@@ -45,7 +54,6 @@ else(OCE_FOUND) #look for OpenCASCADE
/usr/local/opt/opencascade/include
/opt/opencascade/include
/opt/opencascade/inc
/opt/local/include/oce
)
FIND_LIBRARY(OCC_LIBRARY TKernel
/usr/lib
......@@ -56,6 +64,7 @@ else(OCE_FOUND) #look for OpenCASCADE
)
if(OCC_LIBRARY)
message(STATUS "OpenCASCADE has been found.")
GET_FILENAME_COMPONENT(OCC_LIBRARY_DIR ${OCC_LIBRARY} PATH)
IF(NOT OCC_INCLUDE_DIR)
FIND_PATH(OCC_INCLUDE_DIR Standard_Version.hxx
......@@ -112,9 +121,11 @@ if(OCC_FOUND)
TKSTEPAttr
TKHLR
TKFeat
TKXCAF
TKXDESTEP
)
if(OCC_VERSION_STRING VERSION_LESS 6.7)
if(OCC_VERSION_STRING VERSION_LESS 6.8)
MESSAGE(SEND_ERROR "OCC version too low")
endif(OCC_VERSION_STRING VERSION_LESS 6.7)
endif(OCC_VERSION_STRING VERSION_LESS 6.8)
message(STATUS "-- Found OCE/OpenCASCADE with OCC version: ${OCC_VERSION_STRING}")
endif(OCC_FOUND)
......@@ -22,9 +22,34 @@ IF(NEKTAR_USE_MESHGEN)
IF (THIRDPARTY_BUILD_OCE)
INCLUDE(ExternalProject)
SET(OCC_LIBRARIES_TMP PTKernel TKernel TKMath TKBRep TKIGES TKSTEP TKSTEPAttr
TKSTEP209 TKSTEPBase TKShapeSchema TKGeomBase TKGeomAlgo TKG3d TKG2d
TKXSBase TKPShape TKTopAlgo TKShHealing)
SET(OCC_LIBRARIES_TMP
TKFillet
TKMesh
TKernel
TKG2d
TKG3d
TKMath
TKIGES
TKSTL
TKShHealing
TKXSBase
TKBool
TKBO
TKBRep
TKTopAlgo
TKGeomAlgo
TKGeomBase
TKOffset
TKPrim
TKSTEP
TKSTEPBase
TKSTEPAttr
TKHLR
TKFeat
TKXCAF
TKXDESTEP
)
FOREACH(OCC_LIB ${OCC_LIBRARIES_TMP})
LIST(APPEND OCC_LIBRARIES ${TPDIST}/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${OCC_LIB}${CMAKE_SHARED_LIBRARY_SUFFIX})
ENDFOREACH()
......@@ -51,8 +76,6 @@ IF(NEKTAR_USE_MESHGEN)
-DOCE_INSTALL_PREFIX:PATH=${TPDIST}
-DOCE_TESTING=OFF
-DOCE_VISUALISATION=OFF
-DOCE_DISABLE_X11=ON
-DOCE_OCAF=OFF
${TPSRC}/oce-0.17
)
......
......@@ -21,3 +21,11 @@ ADD_NEKTAR_TEST(NodalDemo_Prism_Interp_P7)
ADD_NEKTAR_TEST(NodalDemo_Tet_Deriv_P8)
ADD_NEKTAR_TEST(NodalDemo_Tet_Integral_P6)
ADD_NEKTAR_TEST(NodalDemo_Tet_Interp_P7)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXGear)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXOrder1)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXOrder2)
ADD_NEKTAR_TEST(TimeIntegrationDemoIMEXOrder3)
ADD_NEKTAR_TEST(TimeIntegrationDemoDIRKIMEXOrder2)
ADD_NEKTAR_TEST(TimeIntegrationDemoDIRKIMEXOrder3)
ADD_NEKTAR_TEST(TimeIntegrationDemoCNAB)
ADD_NEKTAR_TEST(TimeIntegrationDemoMCNAB)
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 7</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.10135</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 4</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.0616905</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 5</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.061634</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 6</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.315262</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 1</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.224678</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 2</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.0798986</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 3</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.0610636</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<test>
<description>Test for time integration schemes</description>
<executable>TimeIntegrationDemo</executable>
<parameters>--Npoints 100 --Ntimesteps 100 --NTimeIntegrationMethod 8</parameters>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.134315</value>
</metric>
</metrics>
</test>
......@@ -350,7 +350,16 @@ void Interpolator::Interpolate(
// Obtain Element and LocalCoordinate to interpolate
int elmtid = m_expInField[0]->GetExpIndex(Scoords, Lcoords,
NekConstants::kNekZeroTol);
NekConstants::kGeomFactorsTol);
// we use kGeomFactorsTol as tolerance, while StdPhysEvaluate has
// kNekZeroTol hardcoded, so we need to limit Lcoords to not produce
// a ton of warnings
for(int j = 0; j < nInDim; ++j)
{
Lcoords[j] = std::max(Lcoords[j], -1.0);
Lcoords[j] = std::min(Lcoords[j], 1.0);
}
if (elmtid >= 0)
{
......@@ -431,7 +440,16 @@ void Interpolator::Interpolate(
// Obtain Element and LocalCoordinate to interpolate
int elmtid = m_expInField[0]->GetExpIndex(coords, Lcoords,
NekConstants::kNekZeroTol);
NekConstants::kGeomFactorsTol);
// we use kGeomFactorsTol as tolerance, while StdPhysEvaluate has
// kNekZeroTol hardcoded, so we need to limit Lcoords to not produce
// a ton of warnings
for(int j = 0; j < nInDim; ++j)
{
Lcoords[j] = std::max(Lcoords[j], -1.0);
Lcoords[j] = std::min(Lcoords[j], 1.0);
}
if (elmtid >= 0)
{
......
......@@ -30,7 +30,42 @@
namespace Polylib {
/// The following function is used to circumvent/reduce "Subtractive Cancellation"
/// The expression 1/dz is replaced by optinvsub(.,.)
/// Added on 26 April 2017
double optdiff(double xl, double xr)
{
double m_xln, m_xrn;
int m_expn;
int m_digits = static_cast<int>(fabs(floor(log10(DBL_EPSILON)))-1);
if (fabs(xl-xr)<1.e-4){
m_expn = static_cast<int>(floor(log10(fabs(xl-xr))));
m_xln = xl*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the digits overlap part
m_xrn = xr*powl(10.0L,-m_expn)-floor(xl*powl(10.0L,-m_expn)); // substract the common digits overlap part
m_xln = round(m_xln*powl(10.0L,m_digits+m_expn)); // git rid of rubbish
m_xrn = round(m_xrn*powl(10.0L,m_digits+m_expn));
return powl(10.0L,-m_digits)*(m_xln-m_xrn);
}else{
return (xl-xr);
}
}
double laginterp(double z, int j, const double *zj, int np)
{
double temp = 1.0;
for (int i=0; i<np; i++)
{
if (j != i)
{
temp *=optdiff(z,zj[i])/(zj[j]-zj[i]);
}
}
return temp;
}
/// Define whether to use polynomial deflation (1) or tridiagonal solver (0).
......@@ -1354,29 +1389,15 @@ namespace Polylib {
{
double zi, dz, p, pd, h;
double zi, dz;
zi = *(zgj+i);
dz = z - zi;
if (fabs(dz) < EPS) return 1.0;
jacobd (1, &zi, &pd , np, alpha, beta);
jacobfd(1, &z , &p, NULL , np, alpha, beta);
h = p/(pd*dz);
dz = z-zi;
if (fabs(dz) < EPS) return 1.0;
return h;
return laginterp(z, i, zgj, np);
}
......@@ -1431,34 +1452,15 @@ namespace Polylib {
{
double zi, dz, p, pd, h;
double zi, dz;
zi = *(zgrj+i);
dz = z - zi;
dz = z-zi;
if (fabs(dz) < EPS) return 1.0;
jacobfd (1, &zi, &p , NULL, np-1, alpha, beta + 1);
// need to use this routine in caes zi = -1 or 1
jacobd (1, &zi, &pd, np-1, alpha, beta + 1);
h = (1.0 + zi)*pd + p;
jacobfd (1, &z, &p, NULL, np-1, alpha, beta + 1);
h = (1.0 + z )*p/(h*dz);
return h;
return laginterp(z, i, zgrj, np);
}
......@@ -1515,34 +1517,15 @@ namespace Polylib {
{
double zi, dz, p, pd, h;
double zi, dz;
zi = *(zgrj+i);
dz = z - zi;
dz = z-zi;
if (fabs(dz) < EPS) return 1.0;
jacobfd (1, &zi, &p , NULL, np-1, alpha+1, beta );
// need to use this routine in caes z = -1 or 1
jacobd (1, &zi, &pd, np-1, alpha+1, beta );
h = (1.0 - zi)*pd - p;
jacobfd (1, &z, &p, NULL, np-1, alpha+1, beta);
h = (1.0 - z )*p/(h*dz);
return h;
return laginterp(z, i, zgrj, np);
}
......@@ -1598,35 +1581,15 @@ namespace Polylib {
{
double one = 1., two = 2.;
double zi, dz, p, pd, h;
double zi, dz;
zi = *(zglj+i);
dz = z - zi;
dz = z-zi;
if (fabs(dz) < EPS) return 1.0;
jacobfd(1, &zi, &p , NULL, np-2, alpha + one, beta + one);
// need to use this routine in caes z = -1 or 1
jacobd (1, &zi, &pd, np-2, alpha + one, beta + one);
h = (one - zi*zi)*pd - two*zi*p;
jacobfd(1, &z, &p, NULL, np-2, alpha + one, beta + one);
h = (one - z*z)*p/(h*dz);
return h;
return laginterp(z, i, zglj, np);
}
......
......@@ -742,7 +742,7 @@ namespace Nektar
m_A[0] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
m_B[0] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,twothirdth);
m_A[1] = Array<TwoD,NekDouble>(m_numstages,m_numstages,0.0);
m_B[1] = Array<TwoD,NekDouble>(m_numsteps ,m_numstages,0.0);
m_U = Array<TwoD,NekDouble>(m_numstages,m_numsteps,twothirdth);
m_V = Array<TwoD,NekDouble>(m_numsteps ,m_numsteps, 0.0);
......
......@@ -1488,28 +1488,43 @@ namespace Nektar
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
// Loop over elemente and collect forward expansion
int nexp = GetExpSize();
int n, e, offset, phys_offset;
Array<OneD,NekDouble> e_tmp;
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
{
Vmath::Zero(outarray.num_elements(), outarray, 1);
Array<OneD, NekDouble> tracevals(
m_locTraceToTraceMap->GetNFwdLocTracePts());
m_locTraceToTraceMap->FwdLocTracesFromField(inarray,tracevals);
m_locTraceToTraceMap->
InterpLocEdgesToTrace(0,tracevals,outarray);
m_traceMap->UniversalTraceAssemble(outarray);
}
else
{
// Loop over elemente and collect forward expansion
int nexp = GetExpSize();
int n, e, offset, phys_offset;
Array<OneD,NekDouble> e_tmp;
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
ASSERTL1(outarray.num_elements() >= m_trace->GetNpoints(),
"input array is of insufficient length");
ASSERTL1(outarray.num_elements() >= m_trace->GetNpoints(),
"input array is of insufficient length");
// use m_trace tmp space in element to fill values
for(n = 0; n < nexp; ++n)
{
phys_offset = GetPhys_Offset(n);
for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
// use m_trace tmp space in element to fill values
for(n = 0; n < nexp; ++n)
{
offset = m_trace->GetPhys_Offset(
elmtToTrace[n][e]->GetElmtId());
(*m_exp)[n]->GetEdgePhysVals(e, elmtToTrace[n][e],
inarray + phys_offset,
e_tmp = outarray + offset);
phys_offset = GetPhys_Offset(n);
for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
{
offset = m_trace->GetPhys_Offset(
elmtToTrace[n][e]->GetElmtId());
(*m_exp)[n]->GetEdgePhysVals(e, elmtToTrace[n][e],
inarray + phys_offset,
e_tmp = outarray + offset);
}
}
}
}
......
......@@ -162,6 +162,8 @@ public:
* @brief locates a point in the parametric space
*/
virtual NekDouble loct(Array<OneD, NekDouble> xyz) = 0;
virtual NekDouble DistanceTo(Array<OneD, NekDouble> xyz) = 0;
CADOrientation::Orientation GetOrienationWRT(int surf)
{
......
......@@ -32,15 +32,14 @@
// Description: cad object methods.
//
////////////////////////////////////////////////////////////////////////////////
#include "CADSurf.h"
#include "CADCurve.h"
#include <boost/geometry.hpp>
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include "CADSurf.h"
#include "CADCurve.h"
namespace bg = boost::geometry;
typedef bg::model::d2::point_xy<double> point_xy;
......
......@@ -145,13 +145,17 @@ public:
* @brief takes a point from anywhere find the nearest surface point and its
* uv
*/
virtual void ProjectTo(Array<OneD, NekDouble> &tp,
virtual NekDouble ProjectTo(Array<OneD, NekDouble> &tp,
Array<OneD, NekDouble> &uv) = 0;
virtual Array<OneD, NekDouble> BoundingBox() = 0;
/**
* @brief returns curvature at point uv
*/
virtual NekDouble Curvature(Array<OneD, NekDouble> uv) = 0;
virtual bool IsPlanar() = 0;
/**
* @brief query reversed normal
......@@ -160,6 +164,16 @@ public:
{
return m_orientation;
}
void SetName(std::string i)
{
m_name = i;
}
std::string GetName()
{
return m_name;
}
protected:
/// List of bounding edges in loops with orientation.
......@@ -167,6 +181,8 @@ protected:
/// Function which tests the the value of uv used is within the surface
virtual void Test(Array<OneD, NekDouble> uv) = 0;
std::string m_name;
};
typedef std::shared_ptr<CADSurf> CADSurfSharedPtr;
......
......@@ -101,5 +101,12 @@ Array<OneD, NekDouble> CADSystem::GetPeriodicTranslationVector(int first,
return ret;
}
std::string CADSystem::GetSurfaceName(int i)
{
return m_surfs[i]->GetName();
}
}
}
......@@ -93,7 +93,7 @@ public:
m_2d = false;
}
~CADSystem()
virtual ~CADSystem()
{
}
......@@ -168,8 +168,8 @@ public:
*/
CADCurveSharedPtr GetCurve(int i)