72 logger(
DEBUG,
"LevelSetWall::~LevelSetWall finished");
135 return "LevelSetWall";
197 for (
int i=-
N;
i<
N; ++
i) {
198 for (
int j=-N; j<
N; ++j) {
199 for (
int k=-N; k<
N; ++k) {
201 Vec3D max = min + interactionRadius*Vec3D(2,2,2);
202 Vec3D position = min + interactionRadius*Vec3D(1,1,1);
204 wall.
set(normal,position+distance*normal);
205 std::vector<Vec3D> points;
207 if (!points.empty()) {
221 const unsigned nr = 40;
222 const unsigned nz = 40;
224 std::array<Mdouble, nr> s, c;
225 for (
unsigned i = 0;
i < nr; ++
i)
230 std::array<Mdouble, nz> h, r;
231 for (
unsigned j = 0; j < nz; ++j)
237 for (
unsigned j = 0; j < nz; ++j)
239 for (
unsigned i = 0;
i < nr; ++
i)
246 for (
unsigned j = 0; j < nz - 1; ++j)
248 std::vector<double> cell;
249 cell.reserve(2 * nr + 2);
250 for (
unsigned i = 0;
i < nr; ++
i)
252 cell.push_back(
i + j * nr);
253 cell.push_back(
i + (j + 1) * nr);
255 cell.push_back(j * nr);
256 cell.push_back((j + 1) * nr);
264 s <<
"Values of level set function\n";
267 double meshSize = (
radius_ + extraRadius) / n;
268 for (
int k = -n; k <= n; ++k)
270 for (
int j = -n; j <= n; ++j)
272 for (
int i = -n;
i <= n; ++
i)
284 "data=importdata('LevelSet.txt');\n"
285 "N = nthroot(size(data.data,1),3);\n"
286 "x = reshape(data.data(:,1),N,N,N);\n"
287 "y = reshape(data.data(:,2),N,N,N);\n"
288 "z = reshape(data.data(:,3),N,N,N);\n"
289 "d = reshape(data.data(:,4),N,N,N);\n"
290 "nx= reshape(data.data(:,5),N,N,N);\n"
291 "ny= reshape(data.data(:,6),N,N,N);\n"
292 "nz= reshape(data.data(:,7),N,N,N);\n"
294 "surf(x(:,:,mid),y(:,:,mid),d(:,:,mid),'FaceColor','none')\n"
296 "quiver(x(:,:,mid),y(:,:,mid),nx(:,:,mid),ny(:,:,mid))\n"
297 "legend('distance','normal');\n"
298 "if (min(min(d(:,:,mid)))<0 && max(max(d(:,:,mid)))>0)\n"
299 " contourf(x(:,:,mid),y(:,:,mid),d(:,:,mid),[0 0])\n"
300 " legend('distance','normal','distance>=0'); set(legend,'Location','best');\n"
303 "xlabel('x'); ylabel('y'); title('cross-section at z=0');\n"
304 "view(0,0); axis equal;");
305 logger(
INFO,
"Run LevelSet.m to view level set");
312 Vec3D positionScaled = position /
radius_ *
static_cast<double>(
N) +
Vec3D(1, 1, 1) *
static_cast<double>(
N);
314 int i = std::max(std::min((
int) positionScaled.
X, 2 *
N - 1), 0);
315 int j = std::max(std::min((
int) positionScaled.
Y, 2 *
N - 1), 0);
316 int k = std::max(std::min((
int) positionScaled.
Z, 2 *
N - 1), 0);
318 double x = positionScaled.
X -
i;
319 double y = positionScaled.
Y - j;
320 double z = positionScaled.
Z - k;
322 distance =
levelSet_[
i][j][k] * (1 - x) * (1 - y) * (1 - z)
323 +
levelSet_[i + 1][j][k] * x * (1 - y) * (1 - z)
324 +
levelSet_[i][j + 1][k] * (1 - x) * y * (1 - z)
325 +
levelSet_[i + 1][j + 1][k] * x * y * (1 - z)
326 +
levelSet_[
i][j][k + 1] * (1 - x) * (1 - y) * z
327 +
levelSet_[i + 1][j][k + 1] * x * (1 - y) * z
328 +
levelSet_[i][j + 1][k + 1] * (1 - x) * y * z
329 +
levelSet_[i + 1][j + 1][k + 1] * x * y * z;
330 if (distance < interactionRadius)
333 +
levelSet_[i + 1][j][k] * (1 - y) * (1 - z)
335 +
levelSet_[i + 1][j + 1][k] * y * (1 - z)
337 +
levelSet_[i + 1][j][k + 1] * (1 - y) * z
339 +
levelSet_[i + 1][j + 1][k + 1] * y * z;
342 +
levelSet_[i][j + 1][k] * (1 - x) * (1 - z)
343 +
levelSet_[i + 1][j + 1][k] * x * (1 - z)
346 +
levelSet_[i][j + 1][k + 1] * (1 - x) * z
347 +
levelSet_[i + 1][j + 1][k + 1] * x * z;
353 +
levelSet_[i + 1][j][k + 1] * x * (1 - y)
354 +
levelSet_[i][j + 1][k + 1] * (1 - x) * y
355 +
levelSet_[i + 1][j + 1][k + 1] * x * y;
369 for (
int k = -
N; k <=
N; ++k)
371 for (
int j = -
N; j <=
N; ++j)
373 for (
int i = -
N;
i <=
N; ++
i)
384 for (
int k = -
N; k <=
N; ++k)
386 for (
int j = -
N; j <=
N; ++j)
388 for (
int i = -
N;
i <=
N; ++
i)
401 for (
int k = -
N; k <=
N; ++k)
403 for (
int j = -
N; j <=
N; ++j)
405 for (
int i = -
N;
i <=
N; ++
i)
407 std::array<int, 3> sort;
411 std::sort(sort.begin(), sort.end());
417 else if (sort[0] <
N)
436 for (
int k = -
N; k <=
N; ++k)
438 for (
int j = -
N; j <=
N; ++j)
440 for (
int i = -
N;
i <=
N; ++
i)
451 std::string fileName =
"fourSided.txt";
452 std::ifstream is(fileName.c_str(), std::ios::in);
457 "n=10; % 2n+1 is the number of cells in the level set per dimension\n"
458 "w=1; % [-w,w] is the width of the shape in z-direction\n"
459 "% define a few points on the xy-cross section of the shape that will be interpolated\n"
460 "% we define a quarter of the cross section only; thus the first point has differ \n"
461 "% from the last point by a quarter rotation around the origin\n"
462 "x0=[1 1.8 1.8 1.4 1];\n"
463 "y0=[-1 0 1 1.4 1];\n"
465 "%% create curve, and plot\n"
466 "% define a parameter t0, such that the parametrisation is x0(t0), y0(t0)\n"
468 "% define a more fine-grained parameter t1, and interpolate, obtaining a \n"
469 "% finer parametrisation x1(t1), y1(t1)\n"
470 "t1= linspace(1,length(x0));\n"
471 "x1=interpn(t0,x0,t1,'spline');\n"
472 "y1=interpn(t0,y0,t1,'spline');\n"
473 "% complete the shape by repeating it four times, each rotated by 90 deg\n"
474 "% the points (x,y) are on the shape \n"
475 "x=[x1 -y1 -x1 y1];\n"
476 "y=[y1 x1 -y1 -x1];\n"
477 "% plot the resulting shape\n"
479 "plot(x,y,'-',x0,y0,'o')\n"
482 "%% create level set, and plot\n"
483 "% compute the size of the domain needed to contain the shape ...\n"
484 "radius = max([max(abs(x)),max(abs(y)),abs(w)]);\n"
485 "% .. then divide by the number of cells to get the cell size\n"
486 "radius = n/floor(n/radius);\n"
487 "% create xLS,yLS,zLS values of the level set function dist(xLS,yLS,zLS)\n"
488 "LS = linspace(-radius,radius,2*n+1);\n"
489 "[xLS,yLS,zLS]=meshgrid(LS);\n"
490 "% define the level set function by computing the distance of (xLS,yLS,zLS)\n"
491 "% to the nearest point on the shape (x,y,+-w)\n"
492 "dist = zeros(size(xLS));\n"
493 "for i=1:length(xLS(:))\n"
494 " dist(i) = sqrt(min((x-xLS(i)).^2+(y-yLS(i)).^2));\n"
495 " % make it a signed distance\n"
496 " if inpolygon(xLS(i),yLS(i),x,y)\n"
497 " dist(i) = -dist(i);\n"
499 " % compute the signed distance as the max of signed distances in xy and z\n"
500 " dist(i) = max(dist(i),abs(zLS(i))-w);\n"
502 "% rescale the distance\n"
503 "dist = dist/radius;\n"
505 "% plot the i-th cross-section in xy (i=n+1 is in the center)\n"
508 "mesh(xLS(:,:,i),yLS(:,:,i),dist(:,:,i),'FaceColor','none')\n"
510 "contourf(xLS(:,:,i),yLS(:,:,i),dist(:,:,i),[0 0])\n"
515 "title(['z=' num2str(unique(zLS(:,:,i)))])\n"
516 "xlabel('x'), ylabel('y'), zlabel('dist')\n"
518 "% plot the i-th cross-section in yz (i=n+1 is in the center)\n"
520 "mesh(squeeze(yLS(:,i,:)),squeeze(zLS(:,i,:)),squeeze(dist(:,i,:)),'FaceColor','none')\n"
522 "contourf(squeeze(yLS(:,i,:)),squeeze(zLS(:,i,:)),squeeze(dist(:,i,:)),[0 0])\n"
526 "title(['x=' num2str(unique(xLS(:,i,:)))])\n"
527 "xlabel('y'), ylabel('z'), zlabel('dist')\n"
529 "% plot the i-th cross-section in xz (i=n+1 is in the center)\n"
531 "mesh(squeeze(xLS(i,:,:)),squeeze(zLS(i,:,:)),squeeze(dist(i,:,:)),'FaceColor','none')\n"
533 "contourf(squeeze(xLS(i,:,:)),squeeze(zLS(i,:,:)),squeeze(dist(i,:,:)),[0 0])\n"
537 "title(['y=' num2str(unique(yLS(i,:,:)))])\n"
538 "xlabel('x'), ylabel('z'), zlabel('dist')\n"
540 "%% write n and dist values to file (read into Mercury)\n"
541 "fid = fopen('fourSided.txt','w');\n"
542 "fprintf(fid,'%d\\n',n);\n"
543 "fprintf(fid,'%f ',dist);\n"
545 logger(
ERROR,
"readFromFile: file % could not be opened; make sure you run fourSided.m first.", fileName);
552 logger(
ERROR,
"readFromFile: level set size does not match %", n);
556 for (
int k = -
N; k <=
N; ++k)
558 for (
int j = -
N; j <=
N; ++j)
560 for (
int i = -
N;
i <=
N; ++
i)
void writeToFile(int n, double radiusContact) const
LevelSetWall * copy() const override
Wall copy method. It calls the copy constructor of this Wall, useful for polymorphism.
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
~LevelSetWall() override
Default destructor.
double levelSet_[2 *N+1][2 *N+1][2 *N+1]
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
void createVTK(std::vector< Vec3D > &myPoints) const
const std::complex< Mdouble > i
static void addToVTK(const std::vector< Vec3D > &points, VTKContainer &vtk)
Takes the points provided and adds a triangle strip connecting these points to the vtk container...
void rotate(Vec3D &position) const
Applies the rotation to a position.
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
VTKContainer vtkLabFrame_
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
LevelSetWall(Shape s, double radius, ParticleSpecies *sp=nullptr)
bool getDistanceAndNormalLabCoordinates(Vec3D position, Mdouble interactionRadius, Mdouble &distance, Vec3D &normal) const
std::vector< std::vector< double > > triangleStrips
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
This is a class defining walls.
bool getDistanceAndNormal(const BaseParticle &p, Mdouble &distance, Vec3D &normal_return) const override
Returns the distance of the wall to the particle.
std::vector< Vec3D > points
void read(std::istream &is) override
Reads LevelSetWall from a restart file.
void set(Vec3D normal, Vec3D point)
Defines a standard wall, given an outward normal vector s.t. normal*x=normal*point for all x of the w...
std::string getName() const override
Returns the name of the object, in this case the string "LevelSetWall".
This is a class defining walls.
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
void writeVTK(VTKContainer &vtk) const override
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.