Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.

random surface roughness

Please login with a confirmed email address before reporting spam

In modeling micro-mechanical devices, like switches, surface roughness can have a strong influence. For example, when I have a gap which is closed with a switch, the contact can be strongly affected by the roughness of the relevant surfaces.

I've implemented surface roughness in an external structure generator (Sentaurus Structure Editor), but using external structure generation complicates parameterization. It's nice to be able to form structures within COMSOL.

A displacement field is the obvious approach. A surface roughness function can be defined as a superposition of sine functions:

displacement along axis x = sum over kx { sum over ky { sum over kz { sin(kx * x + dphix_kx,ky,kz) sin(ky * y + dphiy_kx,ky,kz) sin(kz * z + dphiz_kx,ky,kz) Ax_kx_ky_kz }

Here A_kx,ky,kz are amplitudes and dphix_kx,ky,kz are random phases. The amplitudes can be chosen from a normal distribution with zero mean and sigma = exp[-( kx²/2 λx² + ky²/2 λy² + kz²/2 λz² )] where λ's are autocorrelation lengths.

This could be implemented with a loop, but it's not clear how to do that. What would be nice is if a random field was available within comsol, specifying the autocorrelation lengths along each of the three principle axes, and the rms amplitude. Then I could use one of these functions for the relevant displacements. For example, I could specify displacements along just an axis x, or along all three axes x, y, and z, or I could specify such a 3-axis displacement and take the projection to the unit normal.

Is there a way to create rough surfaces within Comsol?

thanks,
Dan

11 Replies Last Post 2016年6月29日 GMT-4 16:55
Jeff Hiller COMSOL Employee

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月13日 GMT-4 10:56
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
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

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月16日 GMT-4 13:23
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.


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: http://www.povray.org/documentation/view/3.6.1/369/ . There's a similar function, "bozo" which is used for clouds.

Jeff Hiller COMSOL Employee

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月17日 GMT-4 12:12
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
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


Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月18日 GMT-4 14:43
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;
}
}
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. http://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 #include #include #include 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 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)

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月20日 GMT-4 16:44
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;
}
}
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 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 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)


Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月20日 GMT-4 19:23
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.
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.

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月21日 GMT-4 09:08
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.
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 Sjodin COMSOL Employee

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月21日 GMT-4 10:27
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.


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 [QUOTE] 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. [/QUOTE]


Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月29日 GMT-4 12:42
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
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

Bjorn Sjodin COMSOL Employee

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月29日 GMT-4 16:50
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
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

Bjorn Sjodin COMSOL Employee

Please login with a confirmed email address before reporting spam

Posted: 8 years ago 2016年6月29日 GMT-4 16:55
I'm attaching an example (created with 5.2a).
I'm attaching an example (created with 5.2a).

Note that while COMSOL employees may participate in the discussion forum, COMSOL® software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.