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 ...


8

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 ...


7

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 $...


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 (...


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 ...


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....


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

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 ...


3

In the light of more recent information: OpenFOAM follows the C++11 standard without any exception at the time of writing. Therefore, you can use any C++ containers of this standard within OpenFOAM. OpenFOAM, at least .com version, has been updated very long time ago to entirely align the ISO/IEC 14882:2011 (i.e. C++11) standard. The main reason why the more ...


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 ...


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 ...


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

All "real" implementations of adaptive mesh refinement start from an unstructured coarse mesh because the domains one wants to solve on are just not always rectangles :-) So your CASE 3 is exactly what people do: in your situation, the coarse mesh just consists of two squares. For this very particular case, one can also implement what is called &...


2

What you're describing is also a critical step in the k-nearest-neighbours method. So no need to reinvent the wheel, we can just look how other people have sped up that algorithm. I don't know about any dictionary like structure that returns this directly, but you could use a k-d tree. If properly implemented, you can get the k closest vectors pretty quickly....


2

You need to look up the VTK file format here: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf It's not very difficult, you'd just write a single cell for each node of your quad tree. The results will look like the pictures you see here or here or here -- all use VTK file format to visualize meshes.


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

Consider the following visualization as an example. It visualizes two binary trees: $T_S$ and $T_V$ for the surface mesh of the sphere and volume mesh of the sphere, respectively. At the 0th level, there is only one node in each tree: $S_1^{(0)}$ and $V_1^{(0)}$. The superscript in the brackets denotes the level in the tree and the subscript denotes the ...


1

As with many questions in computation, a lot of this comes down to what you are trying to achieve. For structured meshes on quadrilaterals, the best and most efficient way to deal with this is not to bother and work implicitly, as you've (kind of suggested). If your cell centred variable has index pair (i,j), then you can label faces to the left/right (...


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 memory ...


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.


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