Source

enzo-dev-problem-types / src / enzo / InitializationKernel.C

/***********************************************************************
/
/  KERNELS FOR INITIALIZATION OF DATA
/
/  written by: Matthew Turk, Devin Silvia
/  date:       January 2012
/
/  PURPOSE:
/
************************************************************************/

#ifdef NEW_PROBLEM_TYPES
#include <string.h>
#include <vector>
#include <stdio.h>
#include <math.h>
#include <iostream>
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
#include "Fluxes.h"
#include "GridList.h"
#include "ExternalBoundary.h"
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
#include "ProblemType.h"
#include "EventHooks.h"

InitializationKernelConstant::InitializationKernelConstant(
    FieldContainer *constants, int additive)
{
    this->constants = constants;
    this->additive = additive;
}

void InitializationKernelConstant::Apply(FieldContainer *fields_,
                FLOAT x, FLOAT y, FLOAT z,
                FLOAT dx, FLOAT dy, FLOAT dz)
{
    FieldContainer fields = (*fields_);
    std::map<std::string, float>::iterator it;
    for(it = this->constants->start();
        it != this->constants->end();
        it++)
    {
      if (this->additive == TRUE) {
        fields[(*it).first] += (*it).second;
      } else {
        fields[(*it).first] = (*it).second;
      }
    }
    return;
}

Enzo::Geometry::SphericalRegion::SphericalRegion(
    FLOAT radius, FLOAT *center)
{
    this->radius = radius;
    this->radius2 = radius * radius;
    this->center[0] = center[0];
    this->center[1] = center[1];
    this->center[2] = center[2];
}

int Enzo::Geometry::SphericalRegion::PointContained(
    FLOAT x, FLOAT y, FLOAT z, FLOAT dx, FLOAT dy, FLOAT dz)
{
    /* If the cell center is within the radius, we return TRUE */
    FLOAT temp[3], r2 = 0;
    int i;
    temp[0] = (x - this->center[0]);
    temp[1] = (y - this->center[1]);
    temp[2] = (z - this->center[2]);
    for (i = 0; i < 3; i++)
        r2 += temp[i] * temp[i];
    if (r2 < this->radius2) return TRUE;
    else return FALSE;
}

Enzo::Geometry::CylindricalRegion::CylindricalRegion(
    FLOAT radius, FLOAT length, FLOAT *center, int axis)
{
    this->radius = radius;
    this->radius2 = radius * radius;
    this->length = length;
    this->halflength = length / 2.0
    this->center[0] = center[0];
    this->center[1] = center[1];
    this->center[2] = center[2];
    this->axis = axis
}

int Enzo::Geometry::CylindericalRegion::PointContained(
    FLOAT x, FLOAT y, FLOAT z, FLOAT dx, FLOAT dy, FLOAT dz)
{
    /* If the cell center is within the radial or lengthwise extent,
    we return TRUE */
    FLOAT temp[3], r2 = 0, l=0;
    int i;
    
    temp[0] = (x - this->center[0]);
    temp[1] = (y - this->center[1]);
    temp[2] = (z - this->center[2]);

    l = abs(temp[this->axis])

    if (axis == 0)
        r2 = (temp[1] * temp[1]) + (temp[2] * temp[2]);

    if (axis == 1)
        r2 = (temp[0] * temp[0]) + (temp[2] * temp[2]);

    if (axis == 2)
        r2 = (temp[0] * temp[0]) + (temp[1] * temp[1]);

    if ( (r2 < this->radius2) && (l < this-halflength) ) return TRUE;
    else return FALSE;
}

#endif