From 30f67d8a0b2c15be88796d5eb303aaed10bd3152 Mon Sep 17 00:00:00 2001
From: Ezgi Berberoglu <berberoglu@biomed.ee.ethz.ch>
Date: Wed, 13 Feb 2019 16:09:31 +0100
Subject: [PATCH] New project function
---
project.py | 121 +++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 121 insertions(+)
create mode 100644 project.py
diff --git a/project.py b/project.py
new file mode 100644
index 0000000..62d5184
--- /dev/null
+++ b/project.py
@@ -0,0 +1,121 @@
+#coding=utf8
+
+################################################################################
+### ###
+### Created by Ezgi BerberoÄŸlu, 2018-2019 ###
+### ###
+### Swiss Federal Institute of Technology (ETH), Zurich, Switzerland ###
+### ###
+### And Martin Genet, 2016-2019 ###
+### ###
+### École Polytechnique, Palaiseau, France ###
+### ###
+################################################################################
+
+import dolfin
+import numpy
+
+import myPythonLibrary as mypy
+import myVTKPythonLibrary as myvtk
+
+import dolfin_dic as ddic
+
+################################################################################
+
+def project(
+ mesh,
+ image_filename,
+ image_field_name="displacement",
+ image_quadrature=1):
+
+ dV = dolfin.Measure("dx", domain=mesh)
+ form_compiler_parameters_for_images = {}
+ form_compiler_parameters_for_images["quadrature_degree"] = image_quadrature
+
+ vfe = dolfin.VectorElement(
+ family="Quadrature",
+ cell=mesh.ufl_cell(),
+ degree=image_quadrature,
+ quad_scheme="default")
+
+ vfs = dolfin.VectorFunctionSpace(
+ mesh=mesh,
+ family="Lagrange",
+ degree=1)
+
+ projected_func = dolfin.Function(
+ vfs,
+ name=image_field_name)
+ U = dolfin.TrialFunction(
+ vfs)
+ V = dolfin.TestFunction(
+ vfs)
+
+ source_expr = dolfin.Expression(
+ cppcode='''\
+#include <vtkSmartPointer.h>
+#include <vtkStructuredGridReader.h>
+#include <vtkProbeFilter.h>
+#include <vtkStructuredGrid.h>
+#include <vtkPolyData.h>
+#include <vtkPointData.h>
+
+namespace dolfin
+{
+
+class MyExpr : public Expression
+{
+ vtkSmartPointer<vtkStructuredGrid> sgrid;
+ vtkSmartPointer<vtkPoints> probe_points;
+ vtkSmartPointer<vtkPolyData> probe_polydata;
+ vtkSmartPointer<vtkProbeFilter> probe_filter;
+
+public:
+
+ MyExpr():
+ Expression(3)
+ {
+ sgrid = vtkSmartPointer<vtkStructuredGrid>::New();
+ probe_points = vtkSmartPointer<vtkPoints>::New();
+ probe_polydata = vtkSmartPointer<vtkPolyData>::New();
+ probe_filter = vtkSmartPointer<vtkProbeFilter>::New();
+ }
+
+ void init_image(
+ const char* image_filename)
+ {
+ vtkSmartPointer<vtkStructuredGridReader> reader = vtkSmartPointer<vtkStructuredGridReader>::New();
+ reader->SetFileName(image_filename);
+ reader->Update();
+ sgrid = reader->GetOutput();
+ probe_filter->SetSourceData(sgrid);
+}
+
+ void eval(
+ Array<double>& expr,
+ const Array<double>& X) const
+ {
+ probe_points->SetNumberOfPoints(1);
+ probe_points->SetPoint(0, X.data());
+ probe_polydata->SetPoints(probe_points);
+ probe_filter->SetInputData(probe_polydata);
+ probe_filter->Update();
+ probe_filter->GetOutput()->GetPointData()->GetArray("'''+image_field_name+'''")->GetTuple(0, expr.data());
+ }
+};
+
+}''',
+ element=vfe)
+ source_expr.init_image(
+ image_filename)
+
+ M = dolfin.assemble(
+ dolfin.inner(U,V)*dV)
+
+ N = dolfin.assemble(
+ dolfin.inner(source_expr,V)*dV,
+ form_compiler_parameters=form_compiler_parameters_for_images)
+
+ dolfin.solve(M, projected_func.vector(), N)
+
+ return projected_func
--
GitLab