Hot answers tagged

7

You have a size mismatch issue: A is a count x 2*count matrix, and you are trying to solve Ax=B with B a 2*count x 1 vector. Moreover, if you compile without -DNDEBUG, you should get a nice assertion telling you this is wrong. To resize a matrix or vector: A.resize(count, 2*count); B.resize(count);


7

The formula still applies for vectors - I find it helpful to use index notation. If $(\nabla u) = u_{,i}$, where $i$ implicitly refers to the $i$th entry of the gradient (i.e. $\frac{\partial u}{\partial x_i}$), and we assume Einstein notation (i.e. repeated indices sum over that index), then we can rewrite $\nabla \cdot \nabla {\bf u} = (u_j)_{,i,i} = \...


6

Like WolfgangBangerth, I strongly recommend that you reconsider your motivation (or your supervisor's motivation) for this goal. First, look at Center for Exascale Simulation of Combustion in Turbulence (ExaCT) web site. They have some proxy applications that may be of use to you in testing your mechanisms, and publications by authors such as John Bell, Joe ...


6

A cross product will tell you both the magnitude and the sign of your angular velocity. In general, angular velocity is defined by a vector as $$\vec{\omega} = \dfrac{\vec{r}\times\vec{v}}{|\vec{r}|^2}$$ where $\vec{\omega}$ is the angular velocity, $\vec{r}$ is the vector from the center of rotation to the point under consideration, and $\vec{v}$ is the ...


4

This is a O(m) problem where m is the size of the big array. So this should be fast: int n = v[v.size()-1]; int k=0; VectorXi lengths[n+1], starts[n+1]; for(int i=0; i<v.size(); ++i) { int i0 = i; while(v[i]==k && i<n) ++i; lengths[k] = i-i0; starts[k] = i; ++k; }


4

What you suggest is a multi-year research program. Don't even think about doing this yourself, you will not within several years achieve anything that comes close to the current research. Rather, search the web for codes that are available and adjust them to your needs. A good first place would be the Department of Energy's Computational Combustion research ...


3

You can calculate the stream function yourself, and from it you can draw contours or streamlines of constant $\psi$. Let us assume you are on a two-dimensional incompressible flow $\mathbf{u}=(u,v,0)^T$, then you can find an exact differential $d\psi$ that satisfies the mass conservation equation: $ \frac{\partial u}{\partial x } + \frac{\partial v}{\...


3

My interpretation of this sentence is the same as yours, and I would have implemented the algorithm the same way.


3

The short answer is "It depends". Certainly, you can try to find the Helmholtz decomposition on your sampled data, and find your irrotational and solenoidal components. However, there are certain requirements on your original vector field you started with. In general, it is that the vector field you are trying to decompose has to be sufficiently smooth and ...


2

For plotting, it is easier in my opinion to not use meshgrid if you want to scale the arrows. You have a vector field $(E_X, E_Z)$ and you can simply normalize it like in the code below: clear; clc; p = 1; e = 8.85*10^(-12); x =linspace(-5 , 5, 20); z = linspace(-5 , 5 ,20); for i = 1:length(x) for k = 1:length(z) R=sqrt(x(i)^2 + z(k)^2) ; ...


2

you just need to add output to your min() function, e.g.: [value, ind] = min(v) ind is then the index values of your min locations (the i,j that you are looking for) and value has the corresponding min values


2

You're correct, if the orientational vector is unitary. Otherwise you must calculate the unitary vector in the direction of $\vec{O}_1$ and then perform the projection (just divide by the norm of the orientational vector). The perpendicular component is just the rejection vector of $\Delta \vec{R}$, given by the definition: $\Delta\vec{R}_{\perp} = \Delta\...


2

I will try to answer the part of the question regarding the accuracy of the calculation, which certainly affects all the other things. Integrals of the type ($T_j$ denoting the $j$th triangle in the discretization and $\vec{x}$ some arbitrary observation point) $$ \int\limits_{T_j}\frac{1}{||\vec{x}-\vec{y}||}dy,\quad \vec{y}\in T_j $$ have a singularity ...


2

You can think of each vector as a point in your linear space. As such, we can use a simple quadtree/octree-like algorithm to map your points into boxes, with "nearby" vectors assigned to the same or an adjacent box. With $n$ total vectors the vector-to-box map costs $\mathcal{O} (n\log n)$, and once this is done you can choose $m$ boxes and select the vector ...


1

For anyone interested in this problem, I found the following solution: #include "itensor/all.h" using namespace itensor; int main() { int N = 100; // // Initialize the site degrees of freedom. // auto sites = SpinHalf(N,{"ConserveQNs=",false}); //make a chain of N spin 1/2's //Transverse field Real h = 4.0; // //...


1

Ok I will answer my own question: The problem is in the line solve(L==a,A,bc) which needs to be replaced by solve(a==L,A,bc). The two versions seem to non interchangeable. Doing so will result in the magnetic field B=project(curl(A),V) which is shown in the picture below:


1

Somehow I have to use t-sne, but I really don't know how. Since you have a PhD my answer will be brief, it's quite a lengthy subject. The aim of dimensionality reduction is to preserve as much of the significant structure of the high-dimensional data as possible in the low-dimensional map. For high-dimensional data that lies on or near a low-dimensional, ...


1

If you know how the x-axis relates from one dataset to the next, the best solution I can think of would be to use the command 'interp1' in matlab. I've compiled a small script with output to show you some example usage. Essentially what you want to do is remap your data from each dataset so that they share some common axis, which you can do with 'interp1' in ...


1

After some searching I found a C/C++ function called project() in the article 'Real-Time Fluid Dynamics for Games' by Jos Stam http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf I had to do a little bit of editing, such as increasing iteration variable k (dependent on how divergent your input field is) and employing periodic boundaries. ...


1

While I don't state this to be an optimal implementation or one taking into account the all edge cases, below is a sample implementation I put together for your task that can give you an idea of how one might approach your problem. The Header File #ifndef _container_hpp_ #define _container_hpp_ #include <vector> class container { public: class ...


1

From your comment "The plot emphasises expectedly the connections of node 1", I guess that maybe the idea is showing that node 1 is connected not only to a set of nodes that are "one close to the other", but its interconnections are spread somewhat evenly across the graph. Reverse Cuthill-McKee in theory reorders the nodes so that clusters are mapped into ...


1

If you just want an approximation to the sum of their norms, you can choose a random sample of $m$ of the indices (i.e., a random set $S \subseteq \{1,2,\dots,N\}$ of size $m$) and compute $$f(\textbf{x}) = \sum_{n \in S} \| x_n \|_2.$$ By a standard central limit theorem argument, $f(\textbf{x})$ will be a reasonable approximation to the sum of the norms; ...


1

Is it possible to tell in index notation whether a vector is a row or column vector, or is that supposed to be clear based on context? It seems the answer is actually lurking in your question itself here: if $u \cdot v \equiv u^T v$ then that doesn't leave much wiggle room. $u^T$ will have to be a row vector in order for the resulting product to be a ...


1

I assume that you are comparing multiplication with an assembled PETSc matrix with your hand-coded matrix-free method. The latter may indeed be faster, but this could be because no entries of A need be loaded from memory. A more meaningful comparison might be to a matrix-free operator in PETSc (see Section 3.3 in the manual). Once you have an apples-to-...


1

You should take a look at the Wikipedia page on rotation matrices, specifically the section on forming a rotation matrix from an axis and an angle. In your case, you have a vector $\vec{v}$ that you would like to align with another vector, $\vec{OP}$. The axis of rotation should be the unit vector $\hat{u}$ normal to these two vectors given by $$\vec{n}=\vec{...


1

I would have expected a geometric discretization to place $E$ and $j$ along edges of the mesh (more precisely, $j$ would be fluxes through facets dual to the mesh edges), in which case $\phi$ is represented by scalar values at the vertices of the mesh. In this way, $E$ and $J$ always "live" at the same mesh elements and are collinear, allowing the dot ...


1

There are several issues that you will need to resolve. 1) In your function alpha_m the loop should be for (int m = 0; m < n; m++) { } because the sum goes from 0 to n-1. 2) I'm not sure what you tried in the function alpha_j. But I believe the loop should not end at m = j-2. Try to match your variable names and starting indices with those of the ...


1

Nvidia's cuSPARSE library implements the functions in the standard sparse blas library, but the sparse blas doesn't include a sparse matrix times sparse vector multiply. However, cuSPARSE also has a sparse matrix times sparse matrix multiplication routine that goes beyond the sparse blas.


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