wip

Makefile

@@ -22,9 +22,11 @@
 ALL: libmbfluid2.so libscontact2.so libscontact3.so
 
 CFLAGS=-Wno-unused-function -g -O3 -march=native -mtune=native -fPIC -std=gnu99 -Iscontact -Ihxt -Ifluid
+CXXFLAGS=-Wno-unused-function -g -O3 -march=native -mtune=native -fPIC -Ihxt -I "/usr/include/eigen3/"
 LDFLAGS=-shared -lm
 
 FLUID_C=fluid/fluid_problem.c fluid/mesh.c fluid/mesh_find.c scontact/quadtree.c hxt/hxt_linear_system.c hxt/hxt_linear_system_lu.c hxt/hxt_message.c fluid/mbxml.c fluid/yxml.c fluid/fluid_problem_io.c
+FLUID_CXX=hxt/hxt_linear_system_eigen.cxx
 FLUID_H=fluid/tools.h fluid/fluid_problem.h hxt/hxt_linear_system.h fluid/mesh_find.h fluid/mesh.h scontact/quadtree.h scontact/vector.h
 
 USE_PETSC?=1
@@ -41,8 +43,9 @@ endif
 SCONTACT_C=scontact/quadtree.c scontact/scontact.c fluid/mbxml.c fluid/yxml.c
 SCONTACT_H=scontact/quadtree.h scontact/scontact.h fluid/mbxml.h fluid/yxml.h
 
-libmbfluid2.so : ${FLUID_C} ${FLUID_H}
-	${CC} ${FLUID_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=2
+libmbfluid2.so : ${FLUID_C} ${FLUID_H} ${FLUID_CXX}
+	${CXX} ${FLUID_CXX} -c -o fluid_cxx.o ${CXXFLAGS}
+	${CC} ${FLUID_C} -o $@ ${CFLAGS} fluid_cxx.o ${LDFLAGS} -DDIMENSION=2
 
 libmbfluid3.so : ${FLUID_C} ${FLUID_H}
 	${CC} ${FLUID_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=3

fluid/fluid_problem.c

@@ -1012,9 +1012,9 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho)
   if (!initialized) {
     hxtInitializeLinearSystems(0, NULL);
     initialized = 1;
-#ifdef HXT_HAVE_PETSC
+/*#ifdef HXT_HAVE_PETSC
     hxtPETScInsertOptions("-pc_type lu -ksp_max_it 30", "fluid");
-#endif
+#endif*/
   }
   FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
   problem->n_fluids = n_fluids;
@@ -1284,11 +1284,12 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
   }
   if (problem->linear_system)
     hxtLinearSystemDelete(&problem->linear_system);
-#ifdef HXT_HAVE_PETSC
+/*#ifdef HXT_HAVE_PETSC
     hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
 #else
     hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
-#endif
+#endif*/
+  hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
   fluid_problem_compute_node_volume(problem);
 }
 
@@ -1368,11 +1369,12 @@ void fluid_problem_after_import(FluidProblem *problem)
   fluid_problem_compute_porosity(problem);
   if(problem->linear_system)
     hxtLinearSystemDelete(&problem->linear_system);
-  #ifdef HXT_HAVE_PETSC
+  /*#ifdef HXT_HAVE_PETSC
     hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements,"fluid");
   #else
     hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
-  #endif 
+  #endif*/
+    hxtLinearSystemCreateEigen(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements);
 }
 
 

hxt/hxt_class_macro.h

@@ -21,7 +21,6 @@
  * see <http://www.gnu.org/licenses/>.
  */
 
-
 #define FE_0(WHAT, _)
 #define FE_1(WHAT, X, _) WHAT(X,1) 
 #define FE_2(WHAT, X, ...) WHAT(X,2)FE_1(WHAT, __VA_ARGS__)

hxt/hxt_linear_system.c

@@ -32,8 +32,6 @@
   operation(opdata,ZeroMatrix,);\
   operation(opdata,Size,int*,);\
   operation(opdata,Solve,double*,double*,);\
-  operation(opdata,SetMatrixRowIdentity,int,int,);\
-  operation(opdata,SetMatrixRowFieldCombinaison,int,int,double*,);\
   operation(opdata,SetRhsEntry,double*,int,int,double,);\
   operation(opdata,AddRhsEntry,double*,int,int,double,);\
   operation(opdata,HasConverged,int*,);\
@@ -43,6 +41,8 @@ HXT_DECLARE_INTERFACE(LinearSystem)
 
 #include "hxt_linear_system_lu.h"
 HXT_DECLARE_DERIVED_CLASS(LinearSystem, LinearSystemLU)
+#include "hxt_linear_system_eigen.h"
+HXT_DECLARE_DERIVED_CLASS(LinearSystem, LinearSystemEigen)
 
 #ifdef HXT_HAVE_PETSC
 #include "hxt_linear_system_petsc.h"
@@ -51,6 +51,7 @@ HXT_DECLARE_DERIVED_CLASS(LinearSystem, LinearSystemPETSc)
 
 HXTStatus hxtInitializeLinearSystems(int *argc, char ***argv) {
   HXT_CHECK(hxtLinearSystemLURegister());
+  HXT_CHECK(hxtLinearSystemEigenRegister());
 #ifdef HXT_HAVE_PETSC
   HXT_CHECK(hxtInitializePETSc(argc, argv));
   HXT_CHECK(hxtLinearSystemPETScRegister());
@@ -67,6 +68,15 @@ HXTStatus hxtLinearSystemCreateLU(HXTLinearSystem **sys, int nElement, int nNode
   return HXT_STATUS_OK;
 };
 
+HXTStatus hxtLinearSystemCreateEigen(HXTLinearSystem **sys, int nElement, int nNodesByElement, int nFields, uint32_t *elements) {
+  if(hxtLinearSystemEigenClass==NULL)
+    HXT_ERROR_MSG(HXT_STATUS_FAILED,"linear system eigen class not registred");
+  HXTLinearSystemEigen *syseigen;
+  HXT_CHECK(hxtLinearSystemEigenCreate(&syseigen, nElement, nNodesByElement, nFields, elements));
+  HXT_CHECK(hxtLinearSystemCreate(sys, hxtLinearSystemEigenClass, syseigen));
+  return HXT_STATUS_OK;
+};
+
 #ifdef HXT_HAVE_PETSC
 HXTStatus hxtLinearSystemCreatePETSc(HXTLinearSystem **sys, int nElement, int nNodesByElement, int nFields, uint32_t *elements, const char *petscOptionsPrefix) {
   if(hxtLinearSystemPETScClass==NULL)

hxt/hxt_linear_system.h

@@ -34,8 +34,6 @@ HXTStatus hxtLinearSystemAddToRhs(HXTLinearSystem *lsys, double *rhs, int el0, c
 HXTStatus hxtLinearSystemSize(HXTLinearSystem *lsys, int *size);
 HXTStatus hxtLinearSystemZeroMatrix(HXTLinearSystem *lsys);
 HXTStatus hxtLinearSystemSolve(HXTLinearSystem *lsys, double *rhs, double *solution);
-HXTStatus hxtLinearSystemSetMatrixRowIdentity(HXTLinearSystem *lsys, int node, int field);
-HXTStatus hxtLinearSystemSetMatrixRowFieldCombinaison(HXTLinearSystem *system, int node, int field, double *coeff);
 HXTStatus hxtLinearSystemSetRhsEntry(HXTLinearSystem *lsys, double *rhs, int node, int field, double v);
 HXTStatus hxtLinearSystemAddRhsEntry(HXTLinearSystem *lsys, double *rhs, int node, int field, double v);
 HXTStatus hxtLinearSystemDelete(HXTLinearSystem **pSystem);
@@ -46,6 +44,8 @@ HXTStatus hxtLinearSystemCreateLU(HXTLinearSystem **sys, int nElement, int nNode
 typedef struct HXTLinearSystemLUStruct HXTLinearSystemLU;
 HXTStatus hxtLinearSystemGetLinearSystemLU(HXTLinearSystem *sys, HXTLinearSystemLU **psys);
 
+HXTStatus hxtLinearSystemCreateEigen(HXTLinearSystem **sys, int nElement, int nNodesByElement, int nFields, uint32_t *elements);
+
 #ifdef HXT_HAVE_PETSC
 typedef struct HXTLinearSystemPETScStruct HXTLinearSystemPETSc;
 HXTStatus hxtLinearSystemCreatePETSc(HXTLinearSystem **sys, int nElement, int nNodesByElement, int nFields, uint32_t *elements, const char *petscOptions);

hxt/hxt_linear_system_eigen.cxx

+/*
+ * MigFlow - Copyright (C) <2010-2018>
+ * <Universite catholique de Louvain (UCL), Belgium
+ *  Universite de Montpellier, France>
+ * 	
+ * List of the contributors to the development of MigFlow: see AUTHORS file.
+ * Description and complete License: see LICENSE file.
+ * 	
+ * This program (MigFlow) is free software: 
+ * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
+ * Public License as published by the Free Software Foundation, either version
+ * 3 of the License, or (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program (see COPYING and COPYING.LESSER files).  If not, 
+ * see <http://www.gnu.org/licenses/>.
+ */
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+extern "C" {
+#include "hxt_linear_system_eigen.h"
+}
+
+#include <Eigen/SparseLU>
+using namespace Eigen;
+
+#define CONMAX 60 
+static void connectivityInsert(int *connectivity, int i, int j)
+{
+  for (int k = 0; k < CONMAX; ++k) {
+    int *p = connectivity + CONMAX*i + k;
+    if (*p == -1)
+      *p = j;
+    if (*p == j)
+      return;
+  }
+  printf("ERROR : node %i has more than %i neighbours\n", i, CONMAX);
+}
+
+struct HXTLinearSystemEigenStruct{
+  SparseMatrix<double,ColMajor> *A;
+  double *x;
+  int *nodeMap;
+  uint32_t *elements;
+  int nNodesByElement;
+  int nElements;
+  int nNodes;
+  int nFields;
+};
+
+static void reverseCuthillMckee(HXTLinearSystemEigen *system, int *ordering)
+{
+  int *nodeConnectivity = (int*)malloc(sizeof(int)*system->nNodes*CONMAX);
+  for (int i = 0; i < system->nNodes*CONMAX; ++i) {
+    nodeConnectivity[i] = -1;
+  }
+  for (int i = 0; i < system->nElements; ++i) {
+    uint32_t *el = system->elements + i*system->nNodesByElement;
+    for (int k = 0; k < system->nNodesByElement; ++k){
+      for (int l = 0; l < system->nNodesByElement; ++l){
+        if (k == l) continue;
+        connectivityInsert(nodeConnectivity, el[k], el[l]);
+        connectivityInsert(nodeConnectivity, el[l], el[k]);
+      }
+    }
+  }
+  int *nodeDegree = (int*)malloc(sizeof(int)*system->nNodes);
+  for (int i = 0; i < system->nNodes; ++i) {
+    nodeDegree[i] = 0;
+    for (int j = 0; j < CONMAX; ++j) {
+      if (nodeConnectivity[CONMAX*i+j] == -1)
+        break;
+      nodeDegree[i] += 1;
+    }
+  }
+  int *queue = (int*)malloc(sizeof(int)*system->nNodes);
+  queue[0] = 0;
+  for (int i = 0; i < system->nNodes; ++i){
+    ordering[i] = -1;
+    if ((nodeDegree[queue[0]] == 0 || nodeDegree[queue[0]] > nodeDegree[i]) && nodeDegree[i] != 0){
+      queue[0] = i;
+    }
+  }
+  int stageStart = 0;
+  int stageEnd = 1;
+  int queueEnd = 1;
+  int id = 0;
+  while(stageStart != stageEnd) {
+    for (int i = stageStart; i < stageEnd; ++i) {
+      int c = queue[i];
+      ordering[c] = id++;
+      for(int j = 0; j < nodeDegree[c]; ++j) {
+        int o = nodeConnectivity[c*CONMAX+j];
+        if (ordering[o] == -1) {
+          ordering[o] = -2;
+#if 1
+          queue[queueEnd++] = o;
+#else
+          int k = stageEnd;
+          while(k < queueEnd && nodeDegree[queue[k]] < nodeDegree[o])
+            k++;
+          for (int l = queueEnd; l > k; --l)
+            queue[l] = queue[l-1];
+          queue[k] = o;
+          queueEnd++;
+#endif
+        }
+      }
+    }
+    stageStart = stageEnd;
+    stageEnd = queueEnd;
+  }
+  for (int i = 0; i < system->nNodes; ++i)
+    if(ordering[i] >= 0)
+      ordering[i] = id-1-ordering[i];
+  free(queue);
+  free(nodeDegree);
+  free(nodeConnectivity);
+}
+
+
+HXTStatus hxtLinearSystemEigenCreate(HXTLinearSystemEigen **pSystem, int nElements, int nNodesByElement, int nFields, uint32_t *elements)
+{
+  HXTLinearSystemEigen *system = (HXTLinearSystemEigen*)malloc(sizeof(HXTLinearSystemEigen));
+  *pSystem = system;
+  system->nFields = nFields;
+  system->nNodesByElement = nNodesByElement;
+  system->nElements = nElements;
+  system->elements = elements;
+  int nNodes = 0;
+  for (int i = 0; i < nElements*nNodesByElement; ++i)
+    if(elements[i]+1 > nNodes)
+      nNodes = elements[i]+1;
+  system->nNodes = nNodes;
+  system->nodeMap = (int*)malloc(sizeof(int)*nNodes);
+  reverseCuthillMckee(system, system->nodeMap);
+  int nUsedNodes = 0;
+  for (int i = 0; i < nNodes; ++i)
+    if (system->nodeMap[i] >= 0)
+      nUsedNodes ++;
+  int n =  nUsedNodes*system->nFields;
+  system->A = new SparseMatrix<double>(n,n);
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenAddToMatrix(HXTLinearSystemEigen *system, int el0, int el1, const double *localMatrix){
+  int nn = system->nNodesByElement;
+  uint32_t *e0 = system->elements + el0*nn;
+  uint32_t *e1 = system->elements + el1*nn;
+  int nf = system->nFields;
+  for (int i = 0; i < nn; ++i) {
+    for (int inf = 0; inf < nf; ++inf) {
+      int ii = system->nodeMap[e0[i]]*nf + inf;
+      for (int j = 0; j < nn; ++j) {
+        for (int jnf = 0; jnf < nf; ++jnf) {
+          int jj = system->nodeMap[e1[j]]*nf + jnf;
+          system->A->coeffRef(ii,jj) += localMatrix[(inf*nn+i)*nf*nn+jnf*nn+j];
+        }
+      }
+    }
+  }
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenAddToRhs(HXTLinearSystemEigen *system, double *rhs, int el0, const double *localVector)
+{
+  int nFields = system->nFields;
+  int nn = system->nNodesByElement;
+  uint32_t *e0 = system->elements + el0*nn;
+  for (int i = 0; i < nFields; ++i) {
+    for (int j = 0; j < nn; ++j) {
+      int m = system->nodeMap[e0[j]]*nFields+i;
+      rhs[m] += localVector[i*nn+j];
+    }
+  }
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenZeroMatrix(HXTLinearSystemEigen *system)
+{
+  printf("Zero !!\n");
+  system->A->setZero();
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenSetRhsEntry(HXTLinearSystemEigen *system, double *rhs, int node, int field, double v)
+{
+  if (system->nodeMap[node] < 0) {
+    return HXT_STATUS_OK;
+  }
+  int row = system->nodeMap[node]*system->nFields + field;
+  rhs[row] = v;
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenAddMatrixEntry(HXTLinearSystemEigen *system, int node0, int field0, int node1, int field1, double v)
+{
+  if (system->nodeMap[node0] < 0 || system->nodeMap[node1] < 0)
+    HXT_ERROR_MSG(HXT_STATUS_FAILED, "node %i or %i not in the domain", node0, node1);
+  int row0 = system->nodeMap[node0]*system->nFields + field0;
+  int col1 = system->nodeMap[node1]*system->nFields + field1;
+  system->A->coeffRef(row0,col1) += v;
+  
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenAddRhsEntry(HXTLinearSystemEigen *system, double *rhs, int node, int field, double v)
+{
+  if (system->nodeMap[node] < 0)
+    return HXT_STATUS_OK;
+  int row = system->nodeMap[node]*system->nFields + field;
+  rhs[row] += v;
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenSolve(HXTLinearSystemEigen *system, double *rhs, double *solution){
+  system->A->makeCompressed();
+  SparseLU<SparseMatrix<double,ColMajor>,COLAMDOrdering<int> > solver;
+  solver.analyzePattern(*system->A);
+  solver.factorize(*system->A);
+  Map<VectorXd> b(rhs,system->A->rows());
+  VectorXd x(solver.solve(b));
+  for (int i = 0; i < system->nNodes; ++i){
+    int ii = system->nodeMap[i];
+    for (int j = 0; j < system->nFields; ++j){
+      solution[i*system->nFields+j] = ii < 0 ? 0. : x[ii*system->nFields+j];
+    }
+  }
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenHasConverged(HXTLinearSystemEigen *lsys, int *converged)
+{
+  *converged = 1;
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenDelete(HXTLinearSystemEigen **pSystem)
+{
+  HXTLinearSystemEigen *system = *pSystem;
+  if (system == NULL)
+    return HXT_STATUS_OK;
+  free(system->A);
+  free(system->x);
+  free(system->nodeMap);
+  free(system);
+  *pSystem = NULL;
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenGetRhsNorm(HXTLinearSystemEigen *lsys, double *rhs, double *pNorm)
+{
+  int size = lsys->A->rows();
+  double norm = 0;
+  for (int i = 0; i < size;++i)
+      norm += rhs[i]*rhs[i];
+  *pNorm =  sqrt(norm);
+  return HXT_STATUS_OK;
+}
+
+HXTStatus hxtLinearSystemEigenSize(HXTLinearSystemEigen *lsys, int *size)
+{
+  *size = lsys->A->rows();
+  return HXT_STATUS_OK;
+}

hxt/hxt_linear_system_eigen.h

+/*
+ * MigFlow - Copyright (C) <2010-2018>
+ * <Universite catholique de Louvain (UCL), Belgium
+ *  Universite de Montpellier, France>
+ * 	
+ * List of the contributors to the development of MigFlow: see AUTHORS file.
+ * Description and complete License: see LICENSE file.
+ * 	
+ * This program (MigFlow) is free software: 
+ * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
+ * Public License as published by the Free Software Foundation, either version
+ * 3 of the License, or (at your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ * 
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program (see COPYING and COPYING.LESSER files).  If not, 
+ * see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef HEXTREME_LINEAR_SYSTEM_EIGEN_H
+#define HEXTREME_LINEAR_SYSTEM_EIGEN_H
+#include <stdint.h>
+#include "hxt_message.h"
+
+typedef struct HXTLinearSystemEigenStruct HXTLinearSystemEigen;
+
+HXTStatus hxtLinearSystemEigenAddToMatrix(HXTLinearSystemEigen *lsys, int el0, int el1, const double *localMatrix);
+HXTStatus hxtLinearSystemEigenAddMatrixEntry(HXTLinearSystemEigen *lsys, int node0, int field0, int node1, int field1, double entry);
+HXTStatus hxtLinearSystemEigenAddToRhs(HXTLinearSystemEigen *lsys, double *rhs, int el0, const double *localVector);
+HXTStatus hxtLinearSystemEigenZeroMatrix(HXTLinearSystemEigen *lsys);
+HXTStatus hxtLinearSystemEigenSolve(HXTLinearSystemEigen *lsys, double *rhs, double *solution);
+HXTStatus hxtLinearSystemEigenSetMatrixRowIdentity(HXTLinearSystemEigen *lsys, int node, int field);
+HXTStatus hxtLinearSystemEigenSetMatrixRowFieldCombinaison(HXTLinearSystemEigen *system, int node, int field, double *coeff);
+HXTStatus hxtLinearSystemEigenSetRhsEntry(HXTLinearSystemEigen *lsys, double *rhs, int node, int field, double v);
+HXTStatus hxtLinearSystemEigenAddRhsEntry(HXTLinearSystemEigen *lsys, double *rhs, int node, int field, double v);
+HXTStatus hxtLinearSystemEigenDelete(HXTLinearSystemEigen **pSystem);
+HXTStatus hxtLinearSystemEigenHasConverged(HXTLinearSystemEigen *lsys, int *converged); 
+HXTStatus hxtLinearSystemEigenGetRhsNorm(HXTLinearSystemEigen *lsys, double *rhs, double *norm);
+HXTStatus hxtLinearSystemEigenCreate(HXTLinearSystemEigen **pSystem, int nElements, int nNodesByElement, int nFields, uint32_t *elements);
+HXTStatus hxtLinearSystemEigenSize(HXTLinearSystemEigen *lsys, int *size);
+
+#endif