ePICURE
Python Isogeometric CUrve REconstruction
This is the final project for the course Applied Mathematics: an Introduction to Numerical Analysis and Sientific Computing. It consists of several small sub-projects, that will combine together in an open source software for curve resconstruction, analysis and simulation. We will apply this to micro-swimming problems, but several other applications are possible...
The students should be familiar with
- Interpolation
- Integration
- Direct and iterative solution of linear systems
- Solution of non-linear systems
- Systems of Ordinary Differential Equations (ODEs)
- Finite difference and Galerkin methods
- Git
Objects
At each time frame, an object can be described by
-
m
numbers (maximum 6), maximum three for the location, and maximum three for the rotation, describing the global position (location and rotation) of the object -
a collection of
n
numbers, describing the local shape of the object. Examples are:- a collection of points, describing a curve
- a collection of points, describing a mesh
- a collection local locations and rotations, describing assemblies of rigid bodies
- use your imagination!
The full time dependent evolution of an object (a stroke) is described by a curve in time for each of the numbers above (a path). The time curves are polynomial curves, or NURBS curves. In any case, they can be described by coefficients times basis functions.
Short Instructions for contributions
Each contribution should be done on a feature branch of a forked repository, then a pull request should be issued, which will be merged by me on mainline. Take a look at this page for a very good explanation on possible git workflows: https://www.atlassian.com/git/tutorials/comparing-workflows/. We will be adopting the forking workflow:
-
Fork the official repository (server side cloning) using the fork button on github
-
Clone your fork onto your PC (it remains yours)
git clone https://github.com/your-user-name/ePICURE.git
-
Add the following master remote to your clone::
cd ePICURE git remote add heltai https://github.com/luca-heltai/ePICURE.git
-
Make a feature branch on your PC
git checkout heltai/master -b a_significant_name_for_your_feature
-
commit, add tests, fix bugs, push your changes to the server, etc. When you push, make sure you push to your remote (
origin
)git push origin a_significant_name_for_your_feature
-
when you think your contribution is ready (a contribution is usually one single feature with one test) do the following:
-
rebase interactively (
-i
switch) your changes onto themaster
branch of the remoteheltai
, after updating your local branches:# Get all changes from all remote repositories git fetch --all git rebase -i heltai/master
-
the above command will pop up an editor in which you will be allowed to make squash all your commits into a couple, more significant, commits. Push your changes by forcing the push (you have just changed history, so this will be allowed only if you force your changes)
git push -f origin a_significant_name_for_your_feature
-
-
go to github web site, select the branch you have been working on, and hit the button compare and pull request
-
once I have reviewed your code, I will accept it an put it in the master branch
Components that should be implemented
Each component should have a well documented interface. An example
is given for the VectorSpace
interface, but all other interfaces
should have an equally detailed interface, which should be discussed
in the issues
section of the git-hub repository.
No contribution will be accepted without a small unit test!. I'm
using nose
to do unit testing:
https://nose.readthedocs.org/en/latest/. In particular it is enough
to write a test file in the directory ./tests. If the filename
contains the word Test
or test
surrounded by _
it will be
executed as a unit. If it contains functions whose name contains
test
or Test
(surrounded by _
or at end/begin of the name)
-
An abstract python interface of type
Shape
, used to describe the shape of a computable object. [TBD: interface] -
An abstract python interface of type
Position
, used to describe the position of a computable object. [TBD: interface] -
An abstract python interface of type
VectorSpace
, used to describe one dimensional functions on[a,b]
, as coefficients times basis functions, with access to single basis functions, their derivatives, their support and the splitting of[a,b]
into sub-intervals where each basis is assumed to be smooth. In particular, the following structure should be implemented for each type of basis:-
VectorSpace.n_dofs
: the number of basis functions that span the space (degrees of freedom) -
VectorSpace.n_dofs_per_end_point
: the number of dofs associated with the end points -
VectorSpace.n_cells
: the number of sub-intervals of[a,b]
where each basis is smooth -
VectorSpace.cells
: then_cells+1
splitting points that make thecells
of theVectorSpace
-
VectorSpace.basis(i)
: the ith basis function (a callable function) -
VectorSpace.basis_der(i, d)
: the d-th derivative of the i-th basis function (a callable function) -
VectorSpace.basis_span(i)
: a tuple indicating the start and end indices into thecells
object where the i-th basis function is different from zero -
VectorSpace.cell_span(i)
: an array of indices containing the basis functions which are non zero on the i-thcell
-
VectorSpace.element(c)
: a callable function, representingsum(c[i] * basis[i])
, which exploits the locality of the basis functions
-
-
The implementation of the above for
- Power basis functions
- Continuous Lagrange (Newton-Cotes with and without end points)
basis functions, of degree
n
with1
repetition - Bernstein basis functions of degree
n
with1
repetition - B-spline basis functions with given degree and knot_vector
- NURBS basis functions with given degree and knot_vector
-
A
CompositeVectorSpace
, taking a subdivision of the interval[a,b]
, aVectorSpace
on[0,1]
and constructing a composite vector space which is a rescaled repetition ofVectorSpace
on each subinterval -
An abstract python interface of type
TimeAnalysis
, capable of taking derivatives and integrals of collections of paths, given an actualVectorSpace
object, and a matrix withVectorSpace.n_dofs
rows and an arbitrary number of columns, representing functions of the givenVectorSpace
-
An
Interpolation
utility, taking aVectorSpace
,n_dofs
points, and constructing a solutionc
to the problemsum_over_j(c[j]*b[j](x[i])) = f(x[i])
-
A
LeastSquare
utility, doing the same but withn>n_dofs
points and using least squares -
A
MassMatrix
assembler, computing the sparse matrix$M_{ij} = \int_0^1 v_i(x) v_j(x) dx$ exploiting locality of theVectorSpace
and aQuadrature
formula given at construction time -
A
StiffnessMatrix
assembler, computing the sparse stiffness matrix$K_{ij} = \int_0^1 v'_i(x) v'_j(x) dx$ exploiting locality of theVectorSpace
and aQuadrature
formula given at construction time -
A
ArcLengthReparametrization
, taking a matrix of coefficients representing a curve in aVectorSpace
, and creating a new one via aLeastSquare
fit which satisfiesc'(x) = L
for eachx
-
A
CurveFromCurvature
function, which takes a vector functiong(x)
representing the curvature of a curve parametrized with arc length, and returns the curve satisfyingc(0) = 0
,c'(x) = L
forx in [0,1]
andc''(x) = g(x)
-
An abstract interface
EquationOfMotion
, taking aShape
and returning aPosition
, such that some underlying equations are satisfied [TBD: Interface]. For this interface, we will construct- collection of sphere swimmers
- filament swimmers using resistive force theory
- crawlers
-
An abstract optimization tool, taking an
EquationOfMotion
, computing aCostFunction
associated to it, imposing aConstraintEquation
on it, and optimizing theShape
evolution, in order to obtain aShape
that satisfies theConstraintEquation
and minimizes theCostFunction
-
A Blender interface
-
given two files, containing numpy nd-arrays (as created by
numpy.save
), representing the time evolution of the position and shape, this interface should create an animation capable of visualising the time evolution of the object. We will need interfaces for the following objects- Golestanian Axis Swimmer (
m=1, n=2
) - Golestanian Plane Swimmer (
m=2, n=3
) - Golestanian Space Swimmer (
m=6, n=4
) - Axis crawler (
m=1, n=as many as you wish
) - Plane worm (
m=2, n=as many as you wish
) - Space worm (
m=6, n=as many as you wish
) - etc.
- Golestanian Axis Swimmer (
-
Naming Conventions
-
Classes and types generally are named using uppercase letters to denote word beginnings (e.g.
VectorSpace
) - sometimes called camel case - while functions and variables use lowercase letters and underscores to separate words -
All file names should use lowercase letters and underscores to separate words
-
Functions or members which return or contain the number of something (number of cells, degrees of freedom, etc) should start with
n_*
. Example:VectorSpace.n_dofs
-
Function which set a bit or flag should start with
set_*
; functions which clear bits of flags should be namedclear_*
-
Function and variable names may not consist of only one or two letters, unless the variable is a pure counting index.
-
Exceptions are used for internal parameter checking and for consistency checks through the
assert
mechanism. -
Each class has to have at least 200 pages of documentation ;-)
Defensive Programming
Defensive programming is a term that we use frequently when we talk about writing code while in the mindset that errors will happen. Here, errors can come in two ways: first, I can make a mistake myself while writing a functions; and secondly, someone else can make a mistake while calling my function. In either case, I would like to write my code in such a way that errors are (i) as unlikely as possible, (ii) that the python interpreter can already find some of the mistakes, and (iii) that the remaining mistakes are relatively easy to find, for example because the program aborts. Defensive programming is then a set of strategies that make these goals more likely.
Over time, we have learned a number of techniques to this end, some of which we list here:
- Assert preconditions on parameters: People call functions with wrong or nonsensical parameters, all the time. Say your function expects to work with two lists of the same length:
def compare_lists(left, right):
# do something with left and right for each element
then make sure that you abort your program if the lists are not compatible:
def compare_lists(left, right):
assert len(left) == len(right)
# do something with left and right for each element
Take a look at https://wiki.python.org/moin/UsingAssertionsEffectively for a more detailed discussion.
Notice that you can turn off assertion automatically: if Python is started with the -O option, then assertions will be stripped out and not evaluated, so you are not adding any overhead to your code. You are only making your life easier...
License:
Please see the file ./LICENSE for details