Merge branch 'cardiac-fibre' of localhost:nektar

CMakeLists.txt

@@ -160,6 +160,7 @@ INCLUDE (ThirdPartyFFTW)
 INCLUDE (ThirdPartyArpack)
 INCLUDE (ThirdPartyMPI)
 INCLUDE (ThirdPartyVTK)
+INCLUDE (ThirdPartyQT4)
 
 IF( NEKTAR_USE_MKL )
     INCLUDE (FindMKL)

cmake/ThirdPartyQT4.cmake

+# QT4 
+OPTION(NEKTAR_USE_QT4
+    "Use QT4 for graphical tools and utilities." OFF)
+
+IF( NEKTAR_USE_QT4 )
+    FIND_PACKAGE(Qt4 COMPONENTS QtCore QtGui REQUIRED)
+    INCLUDE(${QT_USE_FILE})
+    ADD_DEFINITIONS(${QT_DEFINITIONS})
+ENDIF( NEKTAR_USE_QT4 )

dist-exclude

@@ -8,6 +8,7 @@ library/CPackRPM.cmake
 library/Demos/StdRegions/ExtraDemos
 library/Demos/MultiRegions/ExtraDemos
 solvers/ADR2DManifoldSolver
+solvers/CardiacEPSolver/Utilities
 solvers/FitzHughNagumoSolver
 solvers/ImageWarpingSolver
 solvers/PulseWaveSolver

solvers/CardiacEPSolver/CMakeLists.txt

@@ -29,4 +29,6 @@ IF( NEKTAR_SOLVER_CARDIAC_EP )
 
     ADD_SOLVER_EXECUTABLE(CardiacEPSolver solvers-extra 
 			${CardiacEPSolverSource})
+
+    ADD_SUBDIRECTORY(Utilities)
 ENDIF( NEKTAR_SOLVER_CARDIAC_EP )

solvers/CardiacEPSolver/EquationSystems/Monodomain.cpp

@@ -101,10 +101,8 @@ namespace Nektar
                 StdRegions::eVarCoeffD22
         };
         std::string varName = "intensity";
-        std::string varCoeffs[3] = {
-                "AnisotropicConductivityX",
-                "AnisotropicConductivityY",
-                "AnisotropicConductivityZ"
+        std::string aniso_var[3] = {
+                "fx", "fy", "fz"
         };
         int nq = m_fields[0]->GetNpoints();
         Array<OneD, NekDouble> vTemp;
@@ -125,23 +123,9 @@ namespace Nektar
             EvaluateFunction(varName, vTemp, "IsotropicConductivity");
             for (int i = 0; i < m_spacedim; ++i)
             {
-                Vmath::Vmul(nq, vTemp, 1, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
-            }
-        }
-
-        // Apply fibre map (range 0 -> 1)
-        if (m_session->DefinesFunction(varCoeffs[0]))
-        {
-            if (m_session->DefinesCmdLineArgument("verbose"))
-            {
-                cout << "Loading Anisotropic Fibre map." << endl;
-            }
-            for (int i = 0; i < m_spacedim; ++i)
-            {
-                ASSERTL0(m_session->DefinesFunction(varCoeffs[i], varName),
-                    "Function '" + varCoeffs[i] + "' not correctly defined.");
-                EvaluateFunction(varName, vTemp, varCoeffs[i]);
-                Vmath::Vmul(nq, vTemp, 1, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
+                Vmath::Vmul(nq, vTemp, 1,
+                                m_vardiff[varCoeffEnum[i]], 1,
+                                m_vardiff[varCoeffEnum[i]], 1);
             }
         }
 
@@ -186,12 +170,33 @@ namespace Nektar
                                     m_vardiff[varCoeffEnum[i]], 1,
                                     m_vardiff[varCoeffEnum[i]], 1);
                 }
+            }
+        }
+
+        // Apply fibre map (range 0 -> 1)
+        if (m_session->DefinesFunction("AnisotropicConductivity"))
+        {
+            if (m_session->DefinesCmdLineArgument("verbose"))
+            {
+                cout << "Loading Anisotropic Fibre map." << endl;
+            }
+            for (int i = 0; i < m_spacedim; ++i)
+            {
+                ASSERTL0(m_session->DefinesFunction("AnisotropicConductivity",
+                                                    aniso_var[i]),
+                         "Function 'AnisotropicConductivity' not correctly "
+                         "defined.");
+                EvaluateFunction(aniso_var[i], vTemp,
+                                 "AnisotropicConductivity");
+                Vmath::Vmul(nq, vTemp, 1,
+                                m_vardiff[varCoeffEnum[i]], 1,
+                                m_vardiff[varCoeffEnum[i]], 1);
 
                 // Transform variable coefficient and write out to file.
                 m_fields[0]->FwdTrans_IterPerExp(m_vardiff[varCoeffEnum[i]],
                                                  m_fields[0]->UpdateCoeffs());
                 std::stringstream filename;
-                filename << varCoeffs[i];
+                filename << "AnisotropicConductivity_" << aniso_var[i];
                 if (m_comm->GetSize() > 1)
                 {
                     filename << "_P" << m_comm->GetRank();

solvers/CardiacEPSolver/Utilities/CMakeLists.txt

+ADD_SUBDIRECTORY(FibreRegister)
+ADD_SUBDIRECTORY(NavXToVtk)
+ADD_SUBDIRECTORY(FibreToNektar)
+

solvers/CardiacEPSolver/Utilities/FibreRegister/CMakeLists.txt

+IF (NEKTAR_USE_VTK AND NEKTAR_USE_QT4)
+    SET(FR_SOURCES FibreRegister.cpp
+                   MainWindow.cpp)
+
+    SET(QT_OBJECT_HEADERS MainWindow.h)
+
+    QT4_WRAP_CPP(QT_OBJECT_HEADERS_MOC ${QT_OBJECT_HEADERS})
+
+
+    ADD_EXECUTABLE(FibreRegister ${FR_SOURCES} ${QT_OBJECT_HEADERS_MOC})
+    TARGET_LINK_LIBRARIES(FibreRegister ${QT_LIBRARIES})
+    TARGET_LINK_LIBRARIES(FibreRegister vtkCommon vtksys QVTK vtkViews vtkWidgets vtkRendering vtkIO)
+
+    INSTALL(TARGETS FibreRegister 
+	    	RUNTIME DESTINATION ${NEKTAR_BIN_DIR} COMPONENT solvers OPTIONAL)
+ENDIF (NEKTAR_USE_VTK AND NEKTAR_USE_QT4)

solvers/CardiacEPSolver/Utilities/FibreRegister/FibreRegister.cpp

+#include <QtGui/QApplication>
+
+#include "MainWindow.h"
+
+int main(int argc, char* argv[]) {
+    QApplication a( argc, argv );
+    MainWindow w;
+
+    w.show();
+    a.connect( &a, SIGNAL( lastWindowClosed() ), &a, SLOT( quit() ) );
+    return a.exec();
+}

solvers/CardiacEPSolver/Utilities/FibreRegister/MainWindow.h

+#ifndef CLASS_MAINWINDOW_H
+#define CLASS_MAINWINDOW_H
+
+#include <vector>
+#include <string>
+#include <sstream>
+
+#include <QtGui>
+#include <QVTKWidget.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindow.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataReader.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkPointData.h>
+#include <vtkActor.h>
+#include <vtkActor2D.h>
+#include <vtkProperty.h>
+#include <vtkSmoothPolyDataFilter.h>
+#include <vtkDepthSortPolyData.h>
+#include <vtkGlyph3D.h>
+#include <vtkSphereSource.h>
+#include <vtkWorldPointPicker.h>
+#include <vtkIdFilter.h>
+#include <vtkSelectVisiblePoints.h>
+#include <vtkLabeledDataMapper.h>
+#include <vtkContourFilter.h>
+#include <vtkGradientFilter.h>
+#include <vtkPolyDataNormals.h>
+
+class vtkEventQtSlotConnect;
+
+
+class MainWindow : public QMainWindow
+{
+    Q_OBJECT
+
+    public:
+        MainWindow( QWidget* parent = 0, Qt::WindowFlags fl = 0 );
+        ~MainWindow();
+
+    private slots:
+        void BrowseSource();
+        void BrowseTarget();
+
+        void ResetData();
+        void Load();
+        void LoadSource();
+        void LoadSourceLandmarks();
+        void LoadHeightPoints();
+        void LoadTarget();
+
+        void HeightValueChanged(int value);
+        void HeightInterpChanged(int value);
+        void HeightLabelsSelect(bool value);
+
+        void Update();
+        
+        void CreateTargetPoint(vtkObject*, unsigned long, void*, void*, vtkCommand*);
+    
+        void CreateSourcePoint(vtkObject*, unsigned long, void*, void*, vtkCommand*);
+
+        void ExportSource();
+        void ExportTargetPoints();
+        void ExportHeightPoints();
+        void UndoLastLandmarkPoint();
+        void UndoLastHeightPoint();
+
+    private:
+        // GUI widgets
+        QWidget* mRootWidget;
+        QGridLayout* mRootGrid;
+        QGroupBox* mSettingsBox;
+        QVBoxLayout* mSettingsGrid;
+        QGroupBox* mFileBox;
+        QGridLayout* mFileGrid;
+        QLineEdit* mFileSourceEditBox;
+        QPushButton* mFileSourceBrowse;
+        QLineEdit* mFileTargetEditBox;
+        QPushButton* mFileTargetBrowse;
+        QPushButton* mFileLoadButton;
+        QPushButton* mFileExportSourceButton;
+
+        QGroupBox* mLandmarkBox;
+        QGridLayout* mLandmarkGrid;
+        QPushButton* mLandmarksLoadButton;
+        QPushButton* mUndoLandmarkButton;
+        QPushButton* mLandmarksExportButton;
+
+        QGroupBox* mHeightBox;
+        QGridLayout* mHeightGrid;
+        QPushButton* mLoadHeightButton;
+        QPushButton* mUndoHeightButton;
+        QLabel* mHeightSliderLabel;
+        QSlider* mHeightSlider;
+        QSlider* mHeightInterpRange;
+        QCheckBox* mHeightShowLabels;
+        QPushButton* mHeightExport;
+
+        QVTKWidget* mSourceVtk;
+        QVTKWidget* mTargetVtk;
+    
+
+        // VTK
+        // Left Surface pipeline
+        vtkPolyData* mSourceData;
+        vtkLookupTable* mSourceLookupTable;
+        vtkSmoothPolyDataFilter* mSourceFilterSmooth;
+        vtkPolyDataNormals* mSourceNormals;
+        vtkDepthSortPolyData* mSourceFilterDepthSort;
+        vtkPolyDataMapper* mSourceMapper;
+        vtkActor* mSourceActor;
+        vtkRenderer* mSourceRenderer;
+
+        // Right Surface pipeline
+        vtkPolyData* mTargetData;
+        vtkSmoothPolyDataFilter* mTargetFilterSmooth;
+        vtkDepthSortPolyData* mTargetFilterDepthSort;
+        vtkPolyDataMapper* mTargetMapper;
+        vtkActor* mTargetActor;
+        vtkRenderer* mTargetRenderer;
+
+        // Source Points pipeline
+        vtkPolyData* mSourcePointsData;
+        vtkPolyData* mSourceHeightPointData;
+        vtkSphereSource* mSourceSphere;
+        vtkGlyph3D* mSourceFilterGlyph;
+        vtkPolyDataMapper* mSourcePointsMapper;
+        vtkActor* mSourcePointsActor;
+
+        // Source Points Labels pipeline
+        vtkIdFilter* mSourcePointsIds;
+        vtkSelectVisiblePoints* mSourcePointsVisible;
+        vtkSelectVisiblePoints* mSourceHeightPointsVisible;
+        vtkLabeledDataMapper* mSourcePointsLabelMapper;
+        vtkLabeledDataMapper* mSourceHeightPointsLabelMapper;
+        vtkActor2D* mSourcePointsLabelActor;
+        vtkActor2D* mSourceHeightPointsLabelActor;
+    
+        vtkGlyph3D* mSourceHeightGlyph;
+        vtkPolyDataMapper* mSourceHeightPointsMapper;
+        vtkActor* mSourceHeightPointsActor;
+
+        // Source height contours
+        vtkContourFilter* mSourceHeightContours;
+        vtkPolyDataMapper* mSourceHeightContourMapper;
+        vtkActor* mSourceHeightContourActor;
+
+        // Source height gradient
+        vtkPolyDataMapper* mSourceHeightVectorMapper;
+        vtkActor* mSourceHeightVectorActor;
+
+        // Target Points pipeline
+        vtkPolyData* mTargetPointsData;
+        vtkSphereSource* mTargetSphere;
+        vtkGlyph3D* mTargetFilterGlyph;
+        vtkPolyDataMapper* mTargetPointsMapper;
+        vtkActor* mTargetPointsActor;
+
+        // Target Points Labels pipeline
+        vtkIdFilter* mTargetPointsIds;
+        vtkSelectVisiblePoints* mTargetPointsVisible;
+        vtkLabeledDataMapper* mTargetPointsLabelMapper;
+        vtkActor2D* mTargetPointsLabelActor;
+
+        vtkEventQtSlotConnect* Connections;
+        vtkEventQtSlotConnect* Connections_s;
+
+        unsigned int newHeightValue;
+
+        void Draw();
+
+        vtkDoubleArray* InterpSurface(vtkPolyData* pPointData, vtkPolyData* pSurface, int pInterpDistance);
+        void ComputeFibreDirection();
+
+        vtkPoints* AddPoint(vtkPoints* array, double* p);
+        vtkDoubleArray* AddScalar(vtkDataArray* array, double v);
+
+
+};
+
+#endif

solvers/CardiacEPSolver/Utilities/FibreToNektar/CMakeLists.txt

+IF (NEKTAR_USE_VTK)
+    SET(FIBRETONEKTAR_SOURCES FibreToNektar.cpp)
+    ADD_EXECUTABLE(FibreToNektar ${FIBRETONEKTAR_SOURCES})
+    SET_COMMON_PROPERTIES(FibreToNektar)
+    TARGET_LINK_LIBRARIES(FibreToNektar
+            ${NEKTAR++_LIBRARIES}
+            ${Boost_THREAD_LIBRARY}
+            ${Boost_IOSTREAMS_LIBRARY}
+            ${Boost_ZLIB_LIBRARY}
+            ${Boost_DATE_TIME_LIBRARY}
+            ${Boost_PROGRAM_OPTIONS_LIBRARY} 
+            optimized ${TINYXML_LIB} debug ${TINYXML_LIB}
+    )
+    TARGET_LINK_LIBRARIES(FibreToNektar ${TINYXML_LIB})
+    TARGET_LINK_LIBRARIES(FibreToNektar vtkCommon vtksys vtkViews vtkWidgets vtkRendering vtkIO)
+    
+    SET_LAPACK_LINK_LIBRARIES(FibreToNektar)
+            
+    IF (NEKTAR_USE_MPI)
+        TARGET_LINK_LIBRARIES(FibreToNektar ${MPI_LIBRARY} ${MPI_EXTRA_LIBRARY})
+    ENDIF (NEKTAR_USE_MPI)
+    
+    INSTALL(TARGETS FibreToNektar 
+		RUNTIME DESTINATION ${NEKTAR_BIN_DIR} COMPONENT solvers OPTIONAL)
+ENDIF (NEKTAR_USE_VTK)

solvers/CardiacEPSolver/Utilities/FibreToNektar/FibreToNektar.cpp

+#include <LibUtilities/BasicUtils/SessionReader.h>
+#include <LibUtilities/Communication/Comm.h>
+#include <SpatialDomains/MeshGraph2D.h>
+#include <MultiRegions/ExpList2D.h>
+
+#include <vtkPolyDataReader.h>
+#include <vtkPolyData.h>
+#include <vtkPointData.h>
+#include <vtkPoints.h>
+#include <vtkCellArray.h>
+#include <vtkCellDataToPointData.h>
+#include <vtkContourFilter.h>
+// Usage: VtkToFld session.xml input.vtk output.fld
+
+#include <boost/unordered_set.hpp>
+
+
+struct Vertex
+{
+    Vertex(double pX, double pY, double pZ, double pU, double pV, double pW, double factor)
+        : x((int)floor(pX*factor)), y((int)floor(pY*factor)), z((int)floor(pZ*factor)), u(pU), v(pV), w(pW) {}
+    int x;
+    int y;
+    int z;
+    double u;
+    double v;
+    double w;
+
+    bool operator=(const Vertex& v)
+    {
+        return (x == v.x && y == v.y && z == v.z);
+    }
+};
+typedef boost::shared_ptr<Vertex> VertexSharedPtr;
+
+struct VertexHash : std::unary_function<VertexSharedPtr, std::size_t>
+{
+    std::size_t operator()(VertexSharedPtr const& p) const
+    {
+        std::size_t seed = 0;
+        boost::hash_combine(seed, p -> x);
+        boost::hash_combine(seed, p -> y);
+        boost::hash_combine(seed, p -> z);
+        return seed;
+    }
+};
+typedef boost::unordered_set<VertexSharedPtr, VertexHash> VertexSet;
+
+bool operator==(const VertexSharedPtr& v1, const VertexSharedPtr& v2)
+{
+    return v1->operator=(*v2);
+}
+
+
+int main(int argc, char* argv[])
+{
+    if (argc != 5)
+    {
+        cout << "Usage: FibreToNektar [session].xml [fibre].vtk [fibre-out].fld [scale]" << endl;
+        exit(-1);
+    }
+
+    MultiRegions::ExpList2DSharedPtr Exp;
+
+    std::vector<std::string> vFilenames;
+    vFilenames.push_back(std::string(argv[1]));
+
+    LibUtilities::SessionReaderSharedPtr vSession
+            = LibUtilities::SessionReader::CreateInstance(2, argv, vFilenames);
+
+    try
+    {
+        const double factor = atof(argv[4]);
+
+        //----------------------------------------------
+        // Read in mesh from input file
+        SpatialDomains::MeshGraphSharedPtr graph2D = MemoryManager<SpatialDomains::MeshGraph2D>::AllocateSharedPtr(vSession);
+        //----------------------------------------------
+
+        //----------------------------------------------
+        // Define Expansion
+        Exp = MemoryManager<MultiRegions::ExpList2D>::
+            AllocateSharedPtr(vSession,graph2D);
+        //----------------------------------------------
+
+        //----------------------------------------------
+        // Set up coordinates of mesh
+        int coordim = Exp->GetCoordim(0);
+        int nq      = Exp->GetNpoints();
+
+        Array<OneD, NekDouble> xc0(nq,0.0);
+        Array<OneD, NekDouble> xc1(nq,0.0);
+        Array<OneD, NekDouble> xc2(nq,0.0);
+
+        switch(coordim)
+        {
+        case 2:
+            Exp->GetCoords(xc0,xc1);
+            break;
+        case 3:
+            Exp->GetCoords(xc0,xc1,xc2);
+            break;
+        default:
+            ASSERTL0(false,"Coordim not valid");
+            break;
+        }
+        //----------------------------------------------
+
+        vtkPolyDataReader *vtkMeshReader = vtkPolyDataReader::New();
+        vtkMeshReader->SetFileName(argv[2]);
+        vtkMeshReader->Update();
+        vtkPolyData *vtkData = vtkMeshReader->GetOutput();
+
+        vtkPoints *vtkPoints = vtkData->GetPoints();
+
+        ASSERTL0(vtkPoints, "ERROR: cannot get points from mesh.");
+
+        VertexSet points;
+        VertexSet::iterator vIter;
+        double p[3];
+        double vel[3];
+        double x, y, z;
+        int coeff_idx;
+        int i,j;
+        Array<OneD, Array<OneD, NekDouble> > fibre(3);
+        for (i = 0; i < 3; ++i)
+        {
+            fibre[i] = Array<OneD, NekDouble>(Exp->GetNcoeffs());
+        }
+        std::vector<std::string> fibrevar(3);
+        fibrevar[0] = "fx";
+        fibrevar[1] = "fy";
+        fibrevar[2] = "fz";
+
+        // Build up an unordered set of vertices from the VTK file. For each
+        // vertex a hashed value of the coordinates is generated to within a
+        // given tolerance.
+        for (i = 0; i < vtkPoints->GetNumberOfPoints(); ++i)
+        {
+            vtkPoints->GetPoint(i,p);
+            vtkData->GetPointData()->GetArray("FibreDirection")->GetTuple(i, vel);
+            boost::shared_ptr<Vertex> v(new Vertex(p[0],p[1],p[2],vel[0],vel[1],vel[2],factor));
+            points.insert(v);
+        }
+
+        // Now process each vertex of each element in the mesh
+        for (i = 0; i < Exp->GetNumElmts(); ++i)
+        {
+            StdRegions::StdExpansionSharedPtr e = Exp->GetExp(i);
+            for (j = 0; j < e->GetNverts(); ++j)
+            {
+                // Get the index of the coefficient corresponding to this vertex
+                coeff_idx = Exp->GetCoeff_Offset(i) + e->GetVertexMap(j);
+
+                // Get the coordinates of the vertex
+                SpatialDomains::VertexComponentSharedPtr vert = e->GetGeom2D()->GetVertex(j);
+                vert->GetCoords(x,y,z);
+
+                // Look up the vertex in the VertexSet
+                boost::shared_ptr<Vertex> v(new Vertex(x,y,z,0.0, 0.0, 0.0,factor));
+                vIter = points.find(v);
+
+                // If not found, maybe the tolerance should be reduced?
+                // If found, record the scalar value from the VTK file in the
+                // corresponding coefficient.
+                if (vIter == points.end())
+                {
+                    cerr << "Vertex " << i << " not found."
+                            << x << ", " << y << ", " << z << endl;
+                    cerr << v->x << ", " << v->y << ", " << v->z << endl;
+                }
+                else
+                {
+                    double mag = sqrt(((*vIter)->u)*((*vIter)->u) + ((*vIter)->v)*((*vIter)->v) + ((*vIter)->w)*((*vIter)->w));
+                    if (mag < 1e-06) mag = 1.0;
+                    fibre[0][coeff_idx] = fabs((*vIter)->u) / mag;
+                    fibre[1][coeff_idx] = fabs((*vIter)->v) / mag;
+                    fibre[2][coeff_idx] = fabs((*vIter)->w) / mag;
+                }
+            }
+        }
+
+        //-----------------------------------------------
+        // Write solution to file
+        string   out(argv[3]);
+        std::vector<SpatialDomains::FieldDefinitionsSharedPtr> FieldDef
+                                                    = Exp->GetFieldDefinitions();
+        std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
+
+        for (j = 0; j < 3; ++j)
+        {
+            for(i = 0; i < FieldDef.size(); ++i)
+            {
+                FieldDef[i]->m_fields.push_back(fibrevar[j]);
+                Exp->AppendFieldData(FieldDef[i], FieldData[i], fibre[j]);
+            }
+        }
+        graph2D->Write(out, FieldDef, FieldData);
+        //-----------------------------------------------
+    }
+    catch (...) {
+        cout << "An error occurred." << endl;
+    }
+}

solvers/CardiacEPSolver/Utilities/NavXToVtk/CMakeLists.txt

+IF (NEKTAR_USE_VTK)
+    SET(NAVX_SOURCES NavXToVtk.cpp)
+    ADD_EXECUTABLE(NavXToVtk ${NAVX_SOURCES})
+    TARGET_LINK_LIBRARIES(NavXToVtk ${TINYXML_LIB})
+    TARGET_LINK_LIBRARIES(NavXToVtk vtkCommon vtksys vtkViews vtkWidgets vtkRendering vtkIO)
+    INSTALL(TARGETS NavXToVtk 
+		RUNTIME DESTINATION ${NEKTAR_BIN_DIR} COMPONENT solvers OPTIONAL)
+ENDIF (NEKTAR_USE_VTK)

solvers/CardiacEPSolver/Utilities/NavXToVtk/NavXToVtk.cpp

+/*
+ * NavXToVtk.cpp
+ *
+ *  Created on: 19 Jan 2013
+ *      Author: cc
+ */
+
+#include <vtkPolyData.h>
+#include <vtkPolyDataWriter.h>
+#include <vtkDoubleArray.h>