From 5be5c7dd7d1c29381889f10477bf770e6377a1ce Mon Sep 17 00:00:00 2001 From: Spencer Sherwin Date: Sun, 9 Mar 2014 17:12:19 +0000 Subject: [PATCH] First compiling version of InterpField but not debugged --- .../FieldConvert/ProcessInterpField.cpp | 304 ++++++++++++++++++ .../FieldConvert/ProcessInterpField.h | 77 +++++ 2 files changed, 381 insertions(+) create mode 100644 utilities/PostProcessing/FieldConvert/ProcessInterpField.cpp create mode 100644 utilities/PostProcessing/FieldConvert/ProcessInterpField.h diff --git a/utilities/PostProcessing/FieldConvert/ProcessInterpField.cpp b/utilities/PostProcessing/FieldConvert/ProcessInterpField.cpp new file mode 100644 index 000000000..6feb2e35c --- /dev/null +++ b/utilities/PostProcessing/FieldConvert/ProcessInterpField.cpp @@ -0,0 +1,304 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// File: ProcessInterpField.cpp +// +// For more information, please see: http://www.nektar.info/ +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: Interpolate one field to another. +// +//////////////////////////////////////////////////////////////////////////////// +#include +#include +using namespace std; + +#include "ProcessInterpField.h" + +#include +#include + +namespace Nektar +{ + namespace Utilities + { + ModuleKey ProcessInterpField::className = + GetModuleFactory().RegisterCreatorFunction( + ModuleKey(eProcessModule, "interpfield"), + ProcessInterpField::create, "Interpolates one field to another, requires fromxml, fromfld to be defined"); + + + ProcessInterpField::ProcessInterpField(FieldSharedPtr f) : ProcessModule(f) + { + + m_config["fromxml"] = ConfigOption(false,"NotSet","Xml file form which to interpolate field"); + + ASSERTL0(m_config["fromxml"].as().compare("NotSet") != 0, + "Need to specify fromxml=file.xml"); + + m_config["fromfld"] = ConfigOption(false,"NotSet","Fld file form which to interpolate field"); + + ASSERTL0(m_config["fromfld"].as().compare("NotSet") != 0, + "Need to specify fromfld=file.fld "); + + m_config["ClampToLowerValue"] = ConfigOption(false,"-10000", + "Lower bound for interpolation value"); + m_config["ClampToUpperValue"] = ConfigOption(false,"10000", + "Upper bound for interpolation value"); + } + + ProcessInterpField::~ProcessInterpField() + { + } + + void ProcessInterpField::Process(po::variables_map &vm) + { + + if(m_f->m_verbose) + { + cout << "Processing interpolation" << endl; + } + + m_fromField = boost::shared_ptr(new Field()); + + std::vector files; + + // set up session file for from field + files.push_back(m_config["fromxml"].as()); + m_fromField->m_session = LibUtilities::SessionReader:: + CreateInstance(0, 0, files); + + // Set up range based on min and max of local parallel partition + SpatialDomains::DomainRangeShPtr rng = MemoryManager::AllocateSharedPtr(); + + int coordim = m_f->m_exp[0]->GetCoordim(0); + int npts = m_f->m_exp[0]->GetTotPoints(); + Array > coords(3); + + + for(int i = 0; i < coordim; ++i) + { + coords[i] = Array(npts); + } + + for(int i = coordim; i < 3; ++i) + { + coords[i] = NullNekDouble1DArray; + } + + m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]); + + switch(coordim) + { + case 3: + rng->doZrange = true; + rng->zmin = Vmath::Vmin(npts,coords[2],1); + rng->zmax = Vmath::Vmax(npts,coords[2],1); + case 2: + rng->doYrange = true; + rng->ymin = Vmath::Vmin(npts,coords[1],1); + rng->ymax = Vmath::Vmax(npts,coords[1],1); + case 1: + rng->doXrange = true; + rng->xmin = Vmath::Vmin(npts,coords[0],1); + rng->xmax = Vmath::Vmax(npts,coords[0],1); + break; + default: + ASSERTL0(false,"too many values specfied in range"); + } + + // setup rng parameters. + m_fromField->m_graph = SpatialDomains::MeshGraph::Read(m_fromField->m_session,rng); + + // Read in local from field partitions + const SpatialDomains::ExpansionMap &expansions = m_fromField->m_graph->GetExpansions(); + + // check for case where no elements are specified on this + // parallel partition + if(!expansions.size()) + { + return; + } + + Array ElementGIDs(expansions.size()); + SpatialDomains::ExpansionMap::const_iterator expIt; + + int i = 0; + for (expIt = expansions.begin(); expIt != expansions.end(); + ++expIt) + { + ElementGIDs[i++] = expIt->second->m_geomShPtr->GetGlobalID(); + } + + m_fromField->m_fld->Import(m_config["fromfld"].as(), + m_fromField->m_fielddef, + m_fromField->m_data, + LibUtilities::NullFieldMetaDataMap, + ElementGIDs); + + int NumHomogeneousDir = m_f->m_fielddef[0]->m_numHomogeneousDir; + + //---------------------------------------------- + // Set up Expansion information to use mode order from field + m_fromField->m_graph->SetExpansions(m_fromField->m_fielddef); + + int nfields = m_fromField->m_fielddef[0]->m_fields.size(); + + m_fromField->m_exp.resize(nfields); + m_fromField->m_exp[0] = m_f->SetUpFirstExpList(NumHomogeneousDir,true); + + m_f->m_exp.resize(nfields); + + // declare auxiliary fields. + for(i = 1; i < nfields; ++i) + { + m_f->m_exp[i] = m_f->AppendExpList(); + m_fromField->m_exp[i] = m_fromField->AppendExpList(); + } + + + // load field into expansion in fromfield. + for(int j = 0; j < nfields; ++j) + { + for (i = 0; i < m_f->m_fielddef.size(); i++) + { + m_fromField->m_exp[j]->ExtractDataToCoeffs(m_fromField->m_fielddef[i], + m_fromField->m_data[i], + m_fromField->m_fielddef[0]->m_fields[j], + m_fromField->m_exp[j]->UpdateCoeffs()); + } + } + + int nq1 = m_f->m_exp[0]->GetTotPoints(); + + Array x1(nq1); + Array y1(nq1); + Array z1(nq1); + + if (coordim == 2) + { + m_f->m_exp[0]->GetCoords(x1, y1); + } + else if (coordim == 3) + { + m_f->m_exp[0]->GetCoords(x1, y1, z1); + } + + cout << "Interpolating [" << flush; + + NekDouble clamp_low = m_config["ClampToLowerValue"].as(); + NekDouble clamp_up = m_config["ClampToUpperValue"].as(); + InterpolateField(m_fromField->m_exp, m_f->m_exp, + x1, y1, z1, clamp_low, clamp_up); + + cout << "]" << endl; + + + // put field into field data for output + std::vector FieldDef + = m_f->m_exp[0]->GetFieldDefinitions(); + std::vector > FieldData(FieldDef.size()); + + for (int j = 0; j < nfields; ++j) + { + m_f->m_exp[j]->FwdTrans(m_f->m_exp[j]->GetPhys(),m_f->m_exp[j]->UpdateCoeffs()); + for (i = 0; i < FieldDef.size(); ++i) + { + FieldDef[i]->m_fields.push_back(m_fromField->m_fielddef[0]->m_fields[j]); + m_f->m_exp[j]->AppendFieldData(FieldDef[i], FieldData[i]); + } + } + + m_f->m_fielddef = FieldDef; + m_f->m_data = FieldData; + } + + void ProcessInterpField::InterpolateField( + vector &field0, + vector &field1, + Array x, + Array y, + Array z, + NekDouble clamp_low, + NekDouble clamp_up) + { + int expdim = field0[0]->GetCoordim(0); + + Array coords(expdim), Lcoords(expdim); + int nq1 = field1[0]->GetTotPoints(); + int elmtid, offset; + int r, f; + static int intpts = 0; + + ASSERTL0(field0.size() == field1.size(), + "Input field dimension must be same as output dimension"); + + for (r = 0; r < nq1; r++) + { + coords[0] = x[r]; + coords[1] = y[r]; + if (expdim == 3) + { + coords[2] = z[r]; + } + + // Obtain Element and LocalCoordinate to interpolate + elmtid = field0[0]->GetExpIndex(coords, Lcoords, 1e-3); + + offset = field0[0]->GetPhys_Offset(field0[0]-> + GetOffset_Elmt_Id(elmtid)); + + for (f = 0; f < field1.size(); ++f) + { + NekDouble value; + value = field0[f]->GetExp(elmtid)-> + StdPhysEvaluate(Lcoords, field0[f]->GetPhys() +offset); + + if ((boost::math::isnan)(value)) + { + ASSERTL0(false, "new value is not a number"); + } + else + { + value = (value > clamp_up)? clamp_up : + ((value < clamp_low)? clamp_low : + value); + + field1[f]->UpdatePhys()[r] = value; + } + } + + if (intpts%1000 == 0) + { + cout <<"." << flush; + } + intpts ++; + } + } + } +} + + diff --git a/utilities/PostProcessing/FieldConvert/ProcessInterpField.h b/utilities/PostProcessing/FieldConvert/ProcessInterpField.h new file mode 100644 index 000000000..93a738fe8 --- /dev/null +++ b/utilities/PostProcessing/FieldConvert/ProcessInterpField.h @@ -0,0 +1,77 @@ +//////////////////////////////////////////////////////////////////////////////// +// +// File: ProcessInterpField.h +// +// For more information, please see: http://www.nektar.info/ +// +// The MIT License +// +// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), +// Department of Aeronautics, Imperial College London (UK), and Scientific +// Computing and Imaging Institute, University of Utah (USA). +// +// License for the specific language governing rights and limitations under +// Permission is hereby granted, free of charge, to any person obtaining a +// copy of this software and associated documentation files (the "Software"), +// to deal in the Software without restriction, including without limitation +// the rights to use, copy, modify, merge, publish, distribute, sublicense, +// and/or sell copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS +// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +// +// Description: Computes vorticity field. +// +//////////////////////////////////////////////////////////////////////////////// + +#ifndef UTILITIES_PREPROCESSING_FIELDCONVERT_PROCESSVORTICITY +#define UTILITIES_PREPROCESSING_FIELDCONVERT_PROCESSVORTICITY + +#include "Module.h" + +namespace Nektar +{ + namespace Utilities + { + /** + * @brief This processing module interpolates one field to another + */ + class ProcessInterpField : public ProcessModule + { + public: + /// Creates an instance of this class + static boost::shared_ptr create(FieldSharedPtr f) { + return MemoryManager::AllocateSharedPtr(f); + } + static ModuleKey className; + + ProcessInterpField(FieldSharedPtr f); + virtual ~ProcessInterpField(); + + /// Write mesh to output file. + virtual void Process(po::variables_map &vm); + + private: + FieldSharedPtr m_fromField; + + void InterpolateField(vector &field0, + vector &field1, + Array x, + Array y, + Array z, + NekDouble clamp_low, + NekDouble clamp_up); + }; + } +} + +#endif -- GitLab