10

I was pondering this a few days ago (also in Python). Personally I don't think that object oriented programming is always a good fit for numerical programming. You can get distracted with designing the classes rather than just solving the equations. I prefer to stay with simple functions, and with numpy you can have your equations vectorised so the number of ...


9

Simple answer: in modern python every data type is a class, so formally there is no difference between the two solutions you proposed. (Please remember to use new-style classes: classic classes are obsolete! See http://docs.python.org/2/reference/datamodel.html#new-style-and-classic-classes) Now the question should be: how do I organize an efficient data ...


9

I'd suggesting googling for bounding volume hierarchies (BSP tree in particular). Given your point cloud, you can find a plane that splits it into two equal subclouds. Then when you need to find the collection of points that are within some radius R of a test point, you can first compare your test point to that plane, and if it's height above it is more than ...


9

The structure you show is a common choice and equivalent to storing cell-face adjacencies in a CSR matrix format, with the boundary ghost cells in a special place. However, note that FV methods can also be formulated to consist entirely or almost-entirely of face traversal where each face is only visited once (reconstruct to face centroid/quadrature point ...


7

As Aron notes in his comment, rarely are matrices stored in nested pointer data structures in practice. Multiple levels of indirection (I'm told) cause considerable performance penalties, and require multiple allocations and frees (more performance issues, plus more potential to screw up and segfault). I've seen nested pointers used to store matrices only a ...


7

There are several data structure for storing data that preserves information about position and proximity; there by allowing fast nearest neighbor(s) determination. In particular R-trees (and specialized forms like R*-trees) and X-trees. Lots of choices that are optimized for slightly different uses. Choosing a R*-tree rather than a naive nearest neighbor ...


7

Reducing Memory for Sparse Matrices One method (that they mention in the first paper you linked, but is worth emphasizing) is the Block Compressed Sparse Row (BCSR) storage format. If your problem creates dense $n \times n$ blocks (common in e.g. FEM with multiple DoFs per node), you modify the typical CSR (CSC) storage scheme to store only a single column (...


7

One of the authors of fenics, A. Logg, have written a very good paper on datastructures of storing meshes. The paper is A. Logg (2009). Efficient Representation of Computational Meshes http://arxiv.org/abs/1205.3081 In fact it's always a tradeoff between storing all the topological informations (nodes around nodes, faces around nodes, etc...) OR having to ...


6

In deal.II, we basically only use vectors. Maps are too slow and scatter data all around memory, so we typically don't use them if the keys are integers and within a given range. For example, for the connectivity between cells, you can do arrays (STL vectors) in which you store neighbor indices and so that neighbor indices $4i\ldots 4i+3$ correspond to cell $...


6

When I was in the same position as you, I was very happy to find Knuth's Stanford Graphbase (SGB). He supplies not only a library for working with graphs, but also some data sets to play with. My favorite one was the 5-letter English words dataset. You can use your programming language of choice to generate an undirected graph where two words are neighbors ...


5

If you don't want to litter your program with this kind of copying loop then maybe C isn't your language. Assuming this is a common operation, C++ provides far nicer ways to do this sort of thing (by making the array a class and, possibly, making the sub-array a view of it). Even Fortran has a nicer syntax for this kind of operation. What this illustrates is ...


5

This is very similar to one of the biggest challenges in the field of molecular dynamics—computing all of the pairwise interactions between nonbonded particles. There, we use cell lists (or neighbor lists) to help us figure out what's nearby; for this application, the cell list is probably the easier algorithm to use: Divide the box into a series of cells....


5

Although I can't speak for NS or MHD, I do find this "componentization" of function spaces to be a useful design principle in CEM (computational electromagnetics), especially for high-order (in p) discretizations. CEM often uses multiple spaces at the same time: grad-conforming functions to represent electric potential, curl-conforming functions to ...


4

You can do this with any balanced binary search tree data structure by additionally maintaining the total weight below each node. To randomly sample, compute a uniform random number between $0$ and the weight of the root node, and traverse down through the tree until you find the leaf whose range contains the random number. Unfortunately, while maintaining ...


4

This is implemented in qhull which is available from scipy (python). If you cannot use these implementations directly for some reason, the explanations of the data structures in the docs might be helpful. http://www.qhull.org/ http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html#scipy.spatial.Delaunay


4

You should definitely check K-D trees and octrees which are the methods of choice for point sets (while BSPs are for general objects, and grids for more or less uniform densities). They can be very compact and fast, minimizing overhead in both memory and computation, and are simple to implement. When your points are more or less uniformly distributed (even ...


4

deal.II uses variants of your option (2). This is of course also essentially the way you store sparse matrices and is very efficient since you only have two arrays (the array that for each cell stores the beginning index of the data that corresponds to this cell, and the array that contiguously stores all the data of all cells). You traverse both arrays ...


3

I know this question is already answered, but here's a similar single-face based looping storage which is implemented in the OpenFOAM C++ library: Each cell has an index (ID) in a cellList. Two lists are defined for all faces: "face internal owner" and "face neighbour". The length of both face lists corresponds to the number of internal faces in the mesh. ...


3

You should probably consider building the Delaunay triangulation (well, its 3D analogue). In 2D, that's a special triangulation of the data points that always contains the nearest neighbor. The same holds in 3D, but with tetrahedra. You can build once and for all the triangulation, and then search for the nearest neighbor directly in the triangulation. I ...


3

If you order the intervals by their starting point (possible with complexity ${\cal O}(N \log N)$), then you can test whether adjacent intervals overlap and merge them as necessary in ${\cal O}(N)$ operations. What you end up with is list of disjoint intervals. To do this in a periodic domain is only marginally more difficult -- do the same as above, and at ...


3

I suggest to read the Quadpack book (Quadpack, Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2) and Pedro Gonnet's PhD thesis "Adaptive Quadrature Re-revisited (available as pdf here). Pedro is a contributor to ...


3

Straight from Henry's finger tips: a comment which goes in the general direction of the OP's question. As Wolfgang pointed out, OpenFOAM predates many developments we're used to, including the STL and all later developments. Now, the OpenFOAM devs face the question of whether to stick with their implementations and strive for consistency; or to refactor ...


2

Assuming C++, if you're willing to store duplicate items, a time efficient approach for the unweighted case is to store them unsorted in a std::vector or std::deque, v. You can then efficiently add new items by inserting them at the end of v. To draw a sample, simply pick a random index i into v and use sample = v[i]; v[i] = v.back(); v.pop_back(); For the ...


2

you can use http://docs.python.org/2/library/heapq.html as a priority queue


2

I use C++ template to parameterize dimension in a fluid simulation project. It is not directly related to your application, but I think some general ideas would be the same. Vectorize operations for vector/array/point/vertices. I believe you already have this implemented. Provide iterators to access topology in a model. In this way you can abstract(hide) ...


1

The dataset is huge. So I would like to be more memory-efficient. The main point to consider is the following: […] most of the cells will be “unactivated”, i.e most of the elements of the matrix will store an empty list (no signal). This sounds as if an array of lists (or variably sized arrays) is just fine. In this case, the crucial factor for ...


1

Since your dimension is only 3 for the moment, I guess a naive approach would be to create separate trees for each subspace. I would recommend you to use Flann (http://www.cs.ubc.ca/research/flann/) as it is very well implemented and reported to be quite fast for the datasets of your concern (http://www.cs.ubc.ca/research/flann/uploads/FLANN/flann_visapp09....


1

I'm not familiar with the work by Jenny and Lunati, so this answer might be wrong. But if I understand correctly, the geometry you want can be generated in a few lines using PyClaw; see this notebook. This figure shows the primary and dual grids: You can also use mapped grids, and ghost cells are automatically generated if you need them: If either of ...


1

Are you looking for a fast implementation, or are you doing this to learn about Python/numPy/priority queues/Huffman codes? There are a number of different implementations out there already: http://code.activestate.com/recipes/576603-huffman-coding-encoderdeconder/ http://rosettacode.org/wiki/Huffman_coding#Python etc.


1

The data for Algorithms by Sedgewick and Wayne is available at http://introcs.cs.princeton.edu/java/data/.


Only top voted, non community-wiki answers of a minimum length are eligible