Discussion Forum

Hello Daniel,
Whenever I think I need a for loop in pre- or post-processing, my mind goes straight to apps (because apps can call methods that in fact include for loops). Have you considered writing an app with a method that performs that for loop for you?
Jeff

Thanks! I'll look this up in the documentation. However, I think this class of functions is sufficiently useful that it would make sense to provide primary support.

For example, in POV-Ray, an open code used for ray-traced graphics, there is a "bump" function used for surface roughness: www.povray.org/documentation/view/3.6.1/369/ . There's a similar function, "bozo" which is used for clouds.

Hello Daniel,
Depending on what physics are modeled, the impact of surface roughness can sometimes be accounted for without incorporating the roughness into the geometry.
You mentioned ray tracing and, indeed, COMSOL's Ray Tracing Module can account for the effect of surface roughness on rays by introducing a random perturbation to the distribution of the angle at which rays are reflected on a rough surface. See attached screenshot, from the solar dish receiver model in the Ray Optics Module's application library in version 5.2a.
This approach is more economical numerically than modifying the geometry to represent the actual irregularities of the surface.
Best,
Jeff

The principal issue here is the transient vibration behavior of a cantilever. The variation will break symmetry, which allows excitation of antisymmetricmodes, but also result in a difference in the effective spring constant. But you're very correct: there's often tricks that can be done, for example in electrostatic coupling to vary the dielectric constant instead of the thickness of a thin insulating layer.

Anyway, I don't know if this is write, but I tried writing a C function to generate a rough surface. I'm sure there's some errors as it's not debugged fully, although I did compile it as a shared object file and test it with an external C file, plotting the result in gnuplot and checking the statistics (the RMS value is supposed to be 1 but I'm getting closer to 0.8; I'm not sure about why this is).

One issue is the function needs to be initialized with some values (λ, X, Y, Z) and presumably these need to be encoded in the string passed to init but I've not coded this nor worked out the details yet.

This is intended for OS X and Linux, so I've not bothered with any Windows silliness.

i.imgur.com/PcBbC2r.png

/*
surface roughness function of x, y, z
unit auto-correlation length and rms amplitude
use a sprase array of coefficients and phases
*/

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>

float frand() {
return (float)rand()/(float)RAND_MAX;
}

// random normal:
float frandnorm() {
return
frand() + frand() + frand() + frand() +
frand() + frand() + frand() + frand() +
frand() + frand() + frand() + frand() - 6.0;
}

static const char *error = NULL;

const int MAXPOINTS = 1000;
int nkx[MAXPOINTS], nky[MAXPOINTS], nkz[MAXPOINTS];
float a[MAXPOINTS], Φx[MAXPOINTS], Φy[MAXPOINTS], Φz[MAXPOINTS];
const float twoπ = 2 * M_PI;

// parameters for the roughness function
float λ = 5e-9; // autocorrelation length (m)
float X = 1e-6; // dimension of Fourier box (m)
float Y = 1e-6; // dimension of Fourier box (m)
float Z = 1e-6; // dimension of Fourier box (m)

int Nkx = 0; // maximum wavenumber value in direction
int Nky = 0; // maximum wavenumber value in direction
int Nkz = 0; // maximum wavenumber value in direction
int Nk = 0; // total number of k-points
int N = 0;

float kx0, ky0, kz0;

int init(const char *str) {
// kx = 2π N / X, etc
// kmax = 3 / λ
// Nmax = kxmax X / 2π

float kmax = 3.0 / λ; // value of k where attenulation becomes too large
float kmax2 = kmax * kmax;

// increament of k along each axis
kx0 = twoπ / X;
ky0 = twoπ / Y;
kz0 = twoπ / Z;

// number of points along each axis
Nkx = (kx0 == 0) ? MAXPOINTS : floor(1 + kmax / kx0);
Nky = (ky0 == 0) ? MAXPOINTS : floor(1 + kmax / ky0);
Nkz = (kz0 == 0) ? MAXPOINTS : floor(1 + kmax / kz0);

// randomly sample k-space
// populate k-values
if (Nkx * Nky * Nkz < MAXPOINTS) {
int n = 0;
for (int i = 0; i < Nkx; i++) {
float kx = i * kx0;
for (int j = 0; j < Nky; j++) {
float ky = j * ky0;
float k2xy = kx * kx + ky * ky;
if (k2xy > kmax2) continue;
for (int k = 0; k < Nkz; k++) {
float kz = k * kz0;
float k2 = k2xy + kz * kz;
if (k2 > kmax2) continue;
nkx[n] = i;
nky[n] = j;
nkz[n] = k;
n ++;
}
}
}
N = n;
}
else {
// k-space too large: partition space into cubes and randomly sample within each cube
int N1d = 1;
while (N1d * N1d * N1d <= MAXPOINTS)
N1d ++;
N1d --;

if (N1d == 0) N1d = 1;
int n = 0;
for (int i = 0; i < N1d; i ++)
for (int j = 0; j < N1d; j ++)
for (int k = 0; k < N1d; k ++) {
float kx = kmax * (i + frand()) / N1d;
float ky = kmax * (j + frand()) / N1d;
float kz = kmax * (k + frand()) / N1d;
float k2 = kx * kx + ky * ky + kz * kz;
if (k2 > kmax2) continue;
nkx[n] = floor(kx / kx0);
nky[n] = floor(ky / ky0);
nkz[n] = floor(kz / kz0);
n ++;
}
N = n;
}

// assign phases & amplitudes
float sum2 = 0;
for (int i = 0; i < N; i ++) {
Φx[i] = twoπ * frand();
Φy[i] = twoπ * frand();
Φz[i] = twoπ * frand();
float λkx = kx0 * nkx[i] * λ;
float λky = ky0 * nky[i] * λ;
float λkz = kz0 * nkz[i] * λ;
float A = exp(- (λkx * λkx + λky * λky + λkz * λkz) / 2);
a[i] = frandnorm() * A;
sum2 += A * A;
}

// normlize coefficients
if (sum2 > 0) {
float normalization = sqrt(8 / sum2);
for (int i = 0; i < N; i ++)
a[i] *= normalization;
}

return 1;
}

const char * getLastError() {
return error;
}

#define FABS(x) (((x)<0) ? -(x) : x)

int eval(const char *func,
int nArgs,
const double **inReal,
const double **inImag,
int blockSize,
double *outReal,
double *outImag) {
int i, j;
if (strcmp("bump", func) == 0) {
if (nArgs != 3) {
error = "Three arguments expected";
return 0;
}
for (i = 0; i < blockSize; i++) {
double x = inReal[0][i];
double y = inReal[1][i];
double z = inReal[2][i];

/* calculate function */
float f = 0;
for ( int ik = 0; ik < N; ik ++ )
f +=
a[ik] *
(nkx[ik] ? sin(nkx[ik] * kx0 * x + Φx[ik]) : M_SQRT1_2) *
(nky[ik] ? sin(nky[ik] * ky0 * y + Φy[ik]) : M_SQRT1_2) *
(nkz[ik] ? sin(nkx[ik] * kz0 * z + Φz[ik]) : M_SQRT1_2);

outReal[i] = (double) f;
}
return 1;
}
else {
error = "Unknown function";
return 0;
}
}

For linux, I needed to get the unicode specifiers out of that code, and additionally to get rid of some of the C++ line for (int i ...). Anyway, here's the revision, which also calculates derivatives but that's not used since I didn't know how to.

For the interface, I specify a deformed geometry with the following specifications:
x: bump(Xg,Yg,Zg) * 1[nm]
y: bump(Xg+500[nm],Yg,Zg) * 1[nm]
z: 0

This puts bumps in the x-y plane but not z. Ideally I'd have defined separately initialized bump functions for x and y, with separate Fourier coefficients, but as short-hand I offset the values for y well in excess of the autocorrelation length so despite the randomized Fourier coefficients being the same the x and y This is cheating a bit, however.

If the amplitude is increased from 1[nm] to 2[nm] the code fails. It's easy to imagine that if the gradient of a deformation component in the same direction <= -1 then there's a singularity. I suspect that's what's occurring. If I want a greater amplitude I would need to increase the autocorrelation length in conjunction with the amplitude increase to attenuate the high spatial frequency components to prevent this.

The attached image shows the resulting "lumpy" surface. The top surfaces are planar, due to a lack of z-axis displacement.

Note: I've since updated this to accept at most 125 k-points instead of 1000 in attempt to speed it. It's still really slow. The interesting aspect is the deformation is a fixed function of Xg, Yg, Zg and no of time and therefore needs to be calculated only once. I suspect it's being calculated for each point each iteration (I can't be sure of this: I should have the function write to output which should be reflected in the log). One way to do this would be to cache the results. If I had access to a node index that would be easier but I'm not sure it's evaluated only at nodes. I'm not sure...

------------------------------
/*
surface roughness function of x, y, z
unit auto-correlation length and rms amplitude
use a sprase array of coefficients and phases
*/

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

float frand() {
return (float)rand()/(float)RAND_MAX;
}

// random normal:
float frandnorm() {
return
frand() + frand() + frand() + frand() +
frand() + frand() + frand() + frand() +
frand() + frand() + frand() + frand() - 6.0;
}

static const char *error = NULL;

#define MAXPOINTS 1000
int nkx[MAXPOINTS], nky[MAXPOINTS], nkz[MAXPOINTS];
float a[MAXPOINTS], phix[MAXPOINTS], phiy[MAXPOINTS], phiz[MAXPOINTS];
const float twopi = 2 * M_PI;

// parameters for the roughness function
float lambda = 5e-9; // autocorrelation length (m)
float X = 1e-6; // dimension of Fourier box (m)
float Y = 1e-6; // dimension of Fourier box (m)
float Z = 1e-6; // dimension of Fourier box (m)

int Nkx = 0; // maximum wavenumber value in direction
int Nky = 0; // maximum wavenumber value in direction
int Nkz = 0; // maximum wavenumber value in direction
int Nk = 0; // total number of k-points
int N = 0;

float kx0, ky0, kz0;

int init(const char *str) {
// kx = 2pi N / X, etc
// kmax = 3 / lambda
// Nmax = kxmax X / 2pi

// srand((int)time(NULL));

float kmax = 3.0 / lambda; // value of k where attenulation becomes too large
float kmax2 = kmax * kmax;

// increament of k along each axis
kx0 = twopi / X;
ky0 = twopi / Y;
kz0 = twopi / Z;

// number of points along each axis
Nkx = (kx0 == 0) ? MAXPOINTS : floor(1 + kmax / kx0);
Nky = (ky0 == 0) ? MAXPOINTS : floor(1 + kmax / ky0);
Nkz = (kz0 == 0) ? MAXPOINTS : floor(1 + kmax / kz0);

// randomly sample k-space
// populate k-values
if (Nkx * Nky * Nkz < MAXPOINTS) {
int n = 0;
int i;
for (i = 0; i < Nkx; i++) {
float kx = i * kx0;
int j;
for (j = 0; j < Nky; j++) {
float ky = j * ky0;
float k2xy = kx * kx + ky * ky;
if (k2xy > kmax2) continue;
int k;
for (k = 0; k < Nkz; k++) {
float kz = k * kz0;
float k2 = k2xy + kz * kz;
if (k2 > kmax2) continue;
nkx[n] = i;
nky[n] = j;
nkz[n] = k;
n ++;
}
}
}
N = n;
}
else {
// k-space too large: partition space into cubes and randomly sample within each cube
int N1d = 1;
while (N1d * N1d * N1d <= MAXPOINTS)
N1d ++;
N1d --;

if (N1d == 0) N1d = 1;
int n = 0;
int i;
for (i = 0; i < N1d; i ++) {
int j;
for (j = 0; j < N1d; j ++) {
int k;
for (k = 0; k < N1d; k ++) {
float kx = kmax * (i + frand()) / N1d;
float ky = kmax * (j + frand()) / N1d;
float kz = kmax * (k + frand()) / N1d;
float k2 = kx * kx + ky * ky + kz * kz;
if (k2 > kmax2) continue;
nkx[n] = floor(kx / kx0);
nky[n] = floor(ky / ky0);
nkz[n] = floor(kz / kz0);
n ++;
}
}
}
N = n;
}

// assign phases & amplitudes
float sum2 = 0;
int i;
for (i = 0; i < N; i ++) {
phix[i] = twopi * frand();
phiy[i] = twopi * frand();
phiz[i] = twopi * frand();
float lambdakx = kx0 * nkx[i] * lambda;
float lambdaky = ky0 * nky[i] * lambda;
float lambdakz = kz0 * nkz[i] * lambda;
float A = exp(- (lambdakx * lambdakx + lambdaky * lambdaky + lambdakz * lambdakz) / 2);
a[i] = frandnorm() * A;
sum2 += A * A;
}

// normlize coefficients
if (sum2 > 0) {
float normalization = sqrt(8 / sum2);
int i;
for (i = 0; i < N; i ++)
a[i] *= normalization;
}

return 1;
}

const char * getLastError() {
return error;
}

#define FABS(x) (((x)<0) ? -(x) : x)

int eval(const char *func,
int nArgs,
const double **inReal,
const double **inImag,
int blockSize,
double *outReal,
double *outImag) {

if (strcmp("bump", func) == 0) {
if (nArgs != 3) {
error = "Three arguments expected";
return 0;
}
int i;
for (i = 0; i < blockSize; i++) {
double x = inReal[0][i];
double y = inReal[1][i];
double z = inReal[2][i];

/* calculate function */
float f = 0;
int ik;
for (ik = 0; ik < N; ik ++ )
f +=
a[ik] *
(nkx[ik] ? sin(nkx[ik] * kx0 * x + phix[ik]) : M_SQRT1_2) *
(nky[ik] ? sin(nky[ik] * ky0 * y + phiy[ik]) : M_SQRT1_2) *
(nkz[ik] ? sin(nkx[ik] * kz0 * z + phiz[ik]) : M_SQRT1_2);

outReal[i] = (double) f;
}
return 1;
}
/* derivatives:
f = sum { sin(kx x + phix) sin(ky y + phiy) sin(kz z + phiz) }
df/dx = sum { kx cos(kx x + phix) sin(ky y + phiy) sin(kz z + phiz) }
df/dy = sum { ky cos(ky y + phiy) sin(kx x + phix) sin(kz z + phiz) }
df/dz = sum { kz cos(kz z + phiz) sin(kx x + phix) sin(ky y + phiy) }
*/
else {
int deriv = 0;
if (strcmp("dbumpdx", func) == 0)
deriv = 1;
else if (strcmp("dbumpdy", func) == 0)
deriv = 2;
else if (strcmp("dbumpdz", func) == 0)
deriv = 3;
else {
error = "unidentified function";
return 0;
}
if (nArgs != 3) {
error = "Three arguments expected";
return 0;
}
int i;
for (i = 0; i < blockSize; i++) {
double x = inReal[0][i];
double y = inReal[1][i];
double z = inReal[2][i];

/* calculate function */
float f = 0;
int ik;
switch(deriv) {
case 1:
for (ik = 0; ik < N; ik ++ )
if (nkx[ik] != 0)
f +=
a[ik] * nkx[ik] * kx0 *
cos(nkx[ik] * kx0 * x + phix[ik]) *
(nky[ik] ? sin(nky[ik] * ky0 * y + phiy[ik]) : M_SQRT1_2) *
(nkz[ik] ? sin(nkx[ik] * kz0 * z + phiz[ik]) : M_SQRT1_2);
case 2:
for (ik = 0; ik < N; ik ++ )
if (nky[ik] != 0)
f +=
a[ik] * nky[ik] * ky0 *
(nkx[ik] ? sin(nkx[ik] * kx0 * x + phix[ik]) : M_SQRT1_2) *
cos(nky[ik] * ky0 * y + phiy[ik]) *
(nkz[ik] ? sin(nkx[ik] * kz0 * z + phiz[ik]) : M_SQRT1_2);
case 3:
for (ik = 0; ik < N; ik ++ )
if (nkz[ik] != 0)
f +=
a[ik] * nkz[ik] * kz0 *
(nkx[ik] ? sin(nkx[ik] * kx0 * x + phix[ik]) : M_SQRT1_2) *
(nky[ik] ? sin(nky[ik] * ky0 * y + phiy[ik]) : M_SQRT1_2) *
cos(nkz[ik] * kz0 * y + phiz[ik]);
}
outReal[i] = (double) f;
}
return 1;
}
}

I did modify the code to produce independent components (sharing the same k-points, but independent Fourier coefficients) for x, y, and z, so no more kludge shifting coordinates. But it's SLOOOOOW. Clearly the coordinate displacements are being evaluated each iteration, which is unfortunately.

So the key thing now is that I am doing a fairly computation intensive (the weighted sum of 125 sines) for each node, but this is being done each iteration, so this is doing the calculation many, many times for each node, whereas ideally I'd cache it on node vertices, but I'm not sure how to do that.

Alternately in the init method I could create a numerical field in geometry coordinates and then that would be referenced like any other numerical field without the need to recalculate it the hard way.

Another approach would be to find a way to call a function to define a field in definitions, so that it's evaluated only once, rather than each iteration as it is when I refer to my function directly.

Hi Daniel,

I don't think you need to write code to do this. Take a look at the attached screenshot. There you can see how the built-in sum function can be used for this type of task. You can also nest sum functions to create double sums.

Best regards,
Bjorn

So the key thing now is that I am doing a fairly computation intensive (the weighted sum of 125 sines) for each node, but this is being done each iteration, so this is doing the calculation many, many times for each node, whereas ideally I'd cache it on node vertices, but I'm not sure how to do that.

Alternately in the init method I could create a numerical field in geometry coordinates and then that would be referenced like any other numerical field without the need to recalculate it the hard way.

Another approach would be to find a way to call a function to define a field in definitions, so that it's evaluated only once, rather than each iteration as it is when I refer to my function directly.

Bjorn:

Thanks for that! Very nice. What I like about your approach is it avoids the singularity issue with deformation fields when they decrease too rapidly along the axis associated with a given coordinate.

The remaining challenge here is randomizing the roughness. My surface involves 375 randomly chosen amplitudes, 1125 randomly chosen phases, and 1125 randomly chosen k-values. But the key is the function be evaluated only once, and the surface not change over time.

Dan

Hi Dan,

Yes, this could be done as well. You then use an expression like:

scale*sum(sum(if((m!=0)||(n!=0),((m^2+n^2)^(-b/2))*g1(m,n)*cos(2*pi*(m*s1+n*s2+u1(m,n))),0),m,-N,N),n,-N,N)

where:

* scale is a number that determines the approximate max amplitude
* b is a number, typically between 1 and 2. This is a "spectral exponent" determining the "ruggedness" and should be interpreted as 1/f^b where f is the spatial frequency in the spectral distribution of oscillations
* g1(m,n) is a random function with 2 arguments (Defined under Global Definitions>Functions>Random, with Normal distribution (mean 0, standard deviation 1)
* u1(m,n) is another random function with 2 arguments with Uniform distribution (mean 0, standard deviation 1)
* N is the number of spatial frequencies to include

Try, for example:
N=40
scale=0.01
b=1.8
s1 varying from 0.1 to 0.9 (this is really the x-coordinate)
s2 varying from 0.1 to 0.9 (this is really the y-coordinate)
Relative tolerance for the Parametric Surface 1.0E-3
Maximum number of knots for the Parametric Surface 100

The factor ((m^2+n^2)^(-b/2))*g1(m,n) represents a statistical distribution with a power spectrum that decays as 1/f^b. This method is not limited to this type of distribution but you can use any function here f(m,n) to represent a statistical distribution throughout space. You can read more in the book:
"The Science of Fractal Images" by Heinz-Otto Peitgen, Dietmar Saupe

The "if" operator is to avoid division by zero for the "DC" term of the expression.

This is an interesting topic and we may blog about it some day.

I hope this helps.

Bjorn

I'm attaching an example (created with 5.2a).