Commit 6ed257e1 authored by Michael Turner's avatar Michael Turner
Browse files

implementation of nedelar mer optimstaion, does not work

parent 6b02ebf1
......@@ -183,7 +183,7 @@ INCLUDE (ThirdPartyVTK)
INCLUDE (ThirdPartyQT4)
INCLUDE (ThirdPartySMV)
INCLUDE (ThirdPartyTriangle)
INCLUDE (ThirdPartyLevmar)
INCLUDE (ThirdPartyASA047)
INCLUDE (ThirdPartyTetGen)
INCLUDE (ThirdPartyCCM)
......
########################################################################
#
# ThirdParty configuration for Nektar++
#
# Levmar
#
########################################################################
IF(NEKTAR_USE_MESH)
SET(BUILD_LEVMAR ON)
OPTION(THIRDPARTY_BUILD_LEVMAR
"Build Levmar library from ThirdParty." ${BUILD_LEVMAR})
IF (THIRDPARTY_BUILD_LEVMAR)
INCLUDE(ExternalProject)
EXTERNALPROJECT_ADD(
levmar-5.2
PREFIX ${TPSRC}
URL http://ae-nektar.ae.ic.ac.uk/~mt4313/levmar.zip
URL_MD5 7139c9790e3ed4cb5fe2d5be6b1e30d5
STAMP_DIR ${TPBUILD}/stamp
DOWNLOAD_DIR ${TPSRC}
SOURCE_DIR ${TPSRC}/levmar-5.2
BINARY_DIR ${TPBUILD}/levmar-5.2
TMP_DIR ${TPBUILD}/levmar-5.2-tmp
INSTALL_DIR ${TPDIST}
CONFIGURE_COMMAND ${CMAKE_COMMAND}
-G ${CMAKE_GENERATOR}
-DCMAKE_C_COMPILER:FILEPATH=${CMAKE_C_COMPILER}
-DCMAKE_CXX_COMPILER:FILEPATH=${CMAKE_CXX_COMPILER}
-DCMAKE_INSTALL_PREFIX:PATH=${TPDIST}
${TPSRC}/levmar-5.2
)
SET(LEVMAR_LIBRARY levmar CACHE FILEPATH
"Levmar library" FORCE)
SET(LEVMAR_INCLUDE_DIR ${TPDIST}/include CACHE FILEPATH
"Levmar include" FORCE)
LINK_DIRECTORIES(${TPDIST}/lib)
IF (WIN32)
MESSAGE(STATUS
"Build Levmar: ${TPDIST}/${LIB_DIR}/${LEVMAR_LIBRARY}.dll")
ELSE ()
MESSAGE(STATUS
"Build Levmar: ${TPDIST}/${LIB_DIR}/lib${LEVMAR_LIBRARY}.a")
ENDIF ()
SET(TRIANGLE_CONFIG_INCLUDE_DIR ${TPINC})
ELSE()
ADD_CUSTOM_TARGET(levmar-5.2 ALL)
MESSAGE(STATUS "Found Levmar: ${LEVMAR_LIBRARY}")
SET(TRIANGLE_CONFIG_INCLUDE_DIR ${LEVMAR_INCLUDE_DIR})
ENDIF (THIRDPARTY_BUILD_LEVMAR)
INCLUDE_DIRECTORIES(SYSTEM ${LEVMAR_INCLUDE_DIR})
ENDIF(NEKTAR_USE_MESH)
......@@ -91,6 +91,18 @@ Array<OneD, NekDouble> CADSurf::locuv(Array<OneD, NekDouble> p)
return uvr;
}
bool CADSurf::IsPlane()
{
if(m_occSurface.GetType() == GeomAbs_Plane)
{
return true;
}
else
{
return false;
}
}
/**
* @brief Get the x,y,z at parametric point u,v.
*
......
......@@ -95,6 +95,8 @@ class CADSurf
Array<OneD, NekDouble> P (Array<OneD, NekDouble> uv);
Array<OneD, NekDouble> locuv(Array<OneD, NekDouble> p);
bool IsPlane();
private:
/// ID of surface.
int m_ID;
......
......@@ -47,6 +47,7 @@
#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <GeomAbs_SurfaceType.hxx>
#include <BRepTools.hxx>
#include <BRep_Tool.hxx>
#include <gp_Trsf.hxx>
......
......@@ -27,7 +27,7 @@ ADD_NEKTAR_LIBRARY(MeshUtils lib ${NEKTAR_LIBRARY_TYPE}
TARGET_LINK_LIBRARIES(MeshUtils LINK_PUBLIC LibUtilities)
TARGET_LINK_LIBRARIES(MeshUtils LINK_PRIVATE ${TRIANGLE_LIBRARY})
TARGET_LINK_LIBRARIES(MeshUtils LINK_PRIVATE ${TETGEN_LIBRARY})
TARGET_LINK_LIBRARIES(MeshUtils LINK_PRIVATE ${LEVMAR_LIBRARY})
TARGET_LINK_LIBRARIES(MeshUtils LINK_PRIVATE ${ASA_LIBRARY})
INSTALL(DIRECTORY ./
DESTINATION ${NEKTAR_INCLUDE_DIR}/MeshUtils
......
......@@ -34,9 +34,7 @@
////////////////////////////////////////////////////////////////////////////////
#include <MeshUtils/SurfaceMeshing.h>
extern "C"{
#include <levmar.h>
}
#include <asa047.hpp>
using namespace std;
namespace Nektar{
......@@ -44,7 +42,7 @@ namespace MeshUtils {
map<int, MeshNodeSharedPtr> GlobalNodes;
LibUtilities::CADSurfSharedPtr GlobalCad;
int sn, en; //start node end node
int sn, en, m; //start node end node
void SurfaceMeshing::Mesh()
{
......@@ -92,44 +90,47 @@ int sn, en; //start node end node
nodeinlinearmesh = Nodes.size();
}
void EnergyEval(double *p, double *x, int m, int n, void *data)
double EnergyEval(double x[])
{
NekDouble dz = 2.0/(m/2);
Array<OneD, NekDouble> loca,locb;
Array<OneD, NekDouble> uv(2);
x[0] = 0.0;
for(int i = 0; i < m; i++)
{
//cout << x[i] << " ";
}
double p = 0.0;
loca = GlobalNodes[sn]->GetLoc();
uv[0] = p[0]; uv[1] = p[1];
uv[0] = x[0]; uv[1] = x[1];
locb = GlobalCad->P(uv);
x[0] += 1.0/dz*sqrt((loca[0]-locb[0])*(loca[0]-locb[0]) +
p += 1.0/dz*sqrt((loca[0]-locb[0])*(loca[0]-locb[0]) +
(loca[1]-locb[1])*(loca[1]-locb[1]) +
(loca[2]-locb[2])*(loca[2]-locb[2]) );
int i;
for(i = 0; i < m/2 - 1; i++)
{
uv[0] = p[i*2+0]; uv[1] = p[i*2+1];
uv[0] = x[i*2+0]; uv[1] = x[i*2+1];
loca = GlobalCad->P(uv);
uv[0] = p[(i+1)*2+0]; uv[1] = p[(i+1)*2+1];
uv[0] = x[(i+1)*2+0]; uv[1] = x[(i+1)*2+1];
locb = GlobalCad->P(uv);
x[0] += 1.0/dz*sqrt((loca[0]-locb[0])*(loca[0]-locb[0]) +
p += 1.0/dz*sqrt((loca[0]-locb[0])*(loca[0]-locb[0]) +
(loca[1]-locb[1])*(loca[1]-locb[1]) +
(loca[2]-locb[2])*(loca[2]-locb[2]) );
}
uv[0] = p[i*2+0]; uv[1] = p[i*2+1];
uv[0] = x[i*2+0]; uv[1] = x[i*2+1];
loca = GlobalCad->P(uv);
locb = GlobalNodes[en]->GetLoc();
x[0] += 1.0/dz*sqrt((loca[0]-locb[0])*(loca[0]-locb[0]) +
p += 1.0/dz*sqrt((loca[0]-locb[0])*(loca[0]-locb[0]) +
(loca[1]-locb[1])*(loca[1]-locb[1]) +
(loca[2]-locb[2])*(loca[2]-locb[2]) );
for(i = 1; i < m/2; i++)
{
x[i] = 0.0;
}
//cout << x[0] << " " << p[0] << " " << p[1] << endl;
//cout << "\t\t" << p << endl;
return p;
}
void SurfaceMeshing::HOSurf()
......@@ -215,28 +216,56 @@ int sn, en; //start node end node
}
GlobalNodes = Nodes; sn = n[0]; en = n[1]; GlobalCad = m_cad->GetSurf(e->GetSurf());
double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
opts[4]= LM_DIFF_DELTA;
double p[honodes.size()*2];
double x[1];
if(s->IsPlane() == true)
{
e->SetHONodes(honodes);
continue; //optimimum points on plane are linear
}
GlobalNodes = Nodes; sn = n[0]; en = n[1]; GlobalCad = m_cad->GetSurf(e->GetSurf());
m = honodes.size()*2;
cout << "starting" << endl;
double *start;
start = new double[honodes.size()*2];
double *xmin;
xmin = new double[honodes.size()*2];
double *step;
step = new double[honodes.size()*2];
Array<OneD, NekDouble> bound = m_cad->GetSurf(e->GetSurf())->GetBounds();
for(int i = 0; i < honodes.size(); i++)
{
Array<OneD, NekDouble> uv = Nodes[honodes[i]]->GetS(e->GetSurf());
p[i*2+0] = uv[0];
p[i*2+1] = uv[1];
start[i*2+0] = uv[0];
start[i*2+1] = uv[1];
step[i*2+0] = (bound[1]-bound[0])/50.0;
step[i*2+1] = (bound[3]-bound[2])/50.0;
}
double ynew = EnergyEval(start);
//cout << ynew << endl;
int icount, ifault, numres;
nelmin(EnergyEval, honodes.size()*2, start, xmin, &ynew, 1E-8, step,
10, 10000, &icount, &numres, &ifault);
cout << ifault << " " << icount << endl;
if(ifault == 0)
{
for(int i = 0; i < honodes.size(); i++)
{
Array<OneD, NekDouble> uv = Nodes[honodes[i]]->GetS(e->GetSurf());
//cout << uv[0] << " " << uv[1] << endl;
//cout << xmin[i*2+0] << " " << xmin[i*2+1] << endl << endl;
}
//cout << ynew << endl;
//cout << "starting" << endl;
EnergyEval(p,x,honodes.size()*2,honodes.size()*2,NULL);
int ret = dlevmar_dif(&EnergyEval,p,x,honodes.size()*2,honodes.size()*2,1000,opts,info,NULL,NULL,NULL);
for(int i = 0; i < honodes.size(); i++)
{
Array<OneD, NekDouble> uv(2);
uv[0] = xmin[i*2+0]; uv[1] = xmin[i*2+1];
Array<OneD, NekDouble> l = m_cad->GetSurf(e->GetSurf())->P(uv);
Nodes[honodes[i]]->Move(l,uv);
}
}
cout << e->GetId() << " " << ret << endl;
e->SetHONodes(honodes);
}
......
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