Homec4science
Diffusion Dynamic loop scheduling using MPI passive-target remote memory access (master) Browse

dls_MPI_RMA/SpinCalculator.c
4672720499a7master

Last Change
Raw File

SpinCalculator.c
View Options

/*********************************************************************************
* Copyright (c) 2017 *
* Ahmed Eleliemy <ahmed.eleliemy@unibas.ch> *
* University of Basel, Switzerland *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify it *
* under the terms of the license (GNU LGPL) which comes with this package. *
*********************************************************************************/
#include "SpinCalculator.h"
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "SpinDataStructures.h"
#define try do{ jmp_buf ex_buf__; if( !setjmp(ex_buf__) ){
#define catch } else {
#define etry } }while(0)
void calcPointsNormal(MeshObject * mesh)
{
int j=0;
for (j = 0; j < mesh->iPointCount; j++)
{
double cadidpointfacesCount=mesh->facesPerPoint[j];
double candidtevertixNormalX=mesh->dptrXNPoints[j]/cadidpointfacesCount;
double candidtevertixNormalY=mesh->dptrYNPoints[j]/cadidpointfacesCount;
double candidtevertixNormalZ=mesh->dptrZNPoints[j]/cadidpointfacesCount;
double ivalue = sqrt(pow(candidtevertixNormalX, 2) + pow(candidtevertixNormalY, 2) + pow(candidtevertixNormalZ, 2));
candidtevertixNormalX/=(ivalue);
candidtevertixNormalY/=(ivalue);
candidtevertixNormalZ/=(ivalue);
mesh->dptrXNPoints[j]=candidtevertixNormalX;
mesh->dptrYNPoints[j]=candidtevertixNormalY;
mesh->dptrZNPoints[j]=candidtevertixNormalZ;
}
}
void calcSpinImages(double * images,MeshObject * mesh, float binsize, int imagewidth,float supportangle, bool smooth,
int totalSpinpoints, int istartindex,int iendindex)
{
int truth=0;
float supportdistance = binsize * imagewidth;
int iFinishedImages;
int i, k, j, l;
smooth=false;
for (i = istartindex; i < iendindex; i++)
{
double * spinimage = malloc(sizeof(double) * imagewidth * imagewidth);
iFinishedImages=i;
memset(spinimage,0,sizeof(double) *imagewidth*imagewidth);
double centerpointX = mesh->dptrXPoints[i];
double centerpointY = mesh->dptrYPoints[i];
double centerpointZ = mesh->dptrZPoints[i];
double minmumbeta = 0;
double centerpointNormalX=mesh->dptrXNPoints[i];
double centerpointNormalY=mesh->dptrYNPoints[i];
double centerpointNormalZ=mesh->dptrZNPoints[i];
for (j = 0; j < mesh->iPointCount; j++)
{
double candidtevertixX = mesh->dptrXPoints[j];
double candidtevertixY = mesh->dptrYPoints[j];
double candidtevertixZ = mesh->dptrZPoints[j];
double candidtevertixNormalX=mesh->dptrXNPoints[j];
double candidtevertixNormalY=mesh->dptrYNPoints[j];
double candidtevertixNormalZ=mesh->dptrZNPoints[j];
double alpha = 0, beta = 0;
double ab = (centerpointNormalX * candidtevertixNormalX) + (centerpointNormalY * candidtevertixNormalY)
+ (centerpointNormalZ * candidtevertixNormalZ);
double temp=(centerpointNormalX * centerpointNormalX + centerpointNormalY * centerpointNormalY + centerpointNormalZ * centerpointNormalZ);
double a = sqrt(temp);
double b = sqrt((candidtevertixNormalX * candidtevertixNormalX + candidtevertixNormalY * candidtevertixNormalY
+ candidtevertixNormalZ * candidtevertixNormalZ));
if(isnan(a)||isnan(b))
{
truth++;
continue;
}
double angle =acos( ab / (a * b));
if (angle > supportangle)
{
truth++;
continue;
}
double distanceX = candidtevertixX - centerpointX;
double distanceY = candidtevertixY - centerpointY;
double distanceZ = candidtevertixZ - centerpointZ;
beta = centerpointNormalX * distanceX
+ centerpointNormalY * distanceY
+ centerpointNormalZ * distanceZ;
double d_square = distanceX * distanceX
+ distanceY * distanceY + distanceZ * distanceZ;
double beta_square = beta * beta;
alpha = sqrt(d_square - beta_square);
beta = beta - imagewidth * binsize / 2;
if (isnan(alpha)|| beta > 0 || beta < -imagewidth * binsize)
{
truth++;
continue;
}
double indexi = floor(alpha / binsize);
double indexj = floor((-beta) / binsize);
if (indexi >= imagewidth)
{
continue;
}
if (indexj >= imagewidth)
{
continue;
}
a = alpha - indexi * binsize;
b = -beta - indexj * binsize;
k=(int) indexi;
l=(int) indexj;
if(k>imagewidth-1||k<0 || l>imagewidth-1 || l<0 )
{
continue;
}
spinimage[k*imagewidth+l] += 1;
if(k==0 || k==imagewidth-1 || l==0 || l==imagewidth-1)
{
continue;
}
if (smooth)
{
spinimage[(k - 1)*imagewidth + (l - 1)] += a * (1 - b);
spinimage[(k + 1)*imagewidth + (l - 1)] += (1 - a) * (1 - b);
spinimage[(k - 1)*imagewidth + (l + 1)] += a - b;
spinimage[(k + 1)*imagewidth + (l + 1)] += b * (1 - a);
}
}
for (k = 0; k < imagewidth; k++)
{
for (l = 0; l < imagewidth; l++)
{
images[(iFinishedImages*imagewidth*imagewidth)+k*imagewidth+ l]
=ceil(spinimage[k*imagewidth+ l]);
}
}
free(spinimage);
}
}
c4science · Help