C++ DataVector::Initialize方法代码示例

void SolveEvolutionMatrix(
	DataMatrix<double> & matM,
	DataMatrix<double> & matB,
	DataVector<double> & vecAlphaR,
	DataVector<double> & vecAlphaI,
	DataVector<double> & vecBeta,
	DataMatrix<double> & matVR
) {
	int iInfo = 0;

	char jobvl = 'N';
	char jobvr = 'V';

	int n = matM.GetRows();
	int lda = matM.GetRows();
	int ldb = matM.GetRows();
	int ldvl = 1;
	int ldvr = matVR.GetRows();

	DataVector<double> vecWork;
	vecWork.Initialize(8 * n);
	int lwork = vecWork.GetRows();

	dggev_(
		&jobvl,
		&jobvr,
		&n,
		&(matM[0][0]),
		&lda,
		&(matB[0][0]),
		&ldb,
		&(vecAlphaR[0]),
		&(vecAlphaI[0]),
		&(vecBeta[0]),
		NULL,
		&ldvl,
		&(matVR[0][0]),
		&ldvr,
		&(vecWork[0]),
		&lwork,
		&iInfo);

	int nCount = 0;
	for (int i = 0; i < n; i++) {
		if (vecBeta[i] != 0.0) {
			//printf("%i %1.5e %1.5e\n", i, vecAlphaR[i] / vecBeta[i], vecAlphaI[i] / vecBeta[i]);
			nCount++;
		}
	}
/*
	for (int i = 0; i < 40; i++) {
		printf("%1.5e %1.5e\n", matVR[11][4*i+2], matVR[12][4*i+2]);
	}
*/
	Announce("%i total eigenvalues found", nCount);
}

void ReadCFTimeDataFromNcFile(
	NcFile * ncfile,
	const std::string & strFilename,
	std::vector<Time> & vecTimes,
	bool fWarnOnMissingCalendar
) {
	// Empty existing Time vector
	vecTimes.clear();

	// Get time dimension
	NcDim * dimTime = ncfile->get_dim("time");
	if (dimTime == NULL) {
		_EXCEPTION1("Dimension \"time\" not found in file \"%s\"",
			strFilename.c_str());
	}

	// Get time variable
	NcVar * varTime = ncfile->get_var("time");
	if (varTime == NULL) {
		_EXCEPTION1("Variable \"time\" not found in file \"%s\"",
			strFilename.c_str());
	}
	if (varTime->num_dims() != 1) {
		_EXCEPTION1("Variable \"time\" has more than one dimension in file \"%s\"",
			strFilename.c_str());
	}
	if (strcmp(varTime->get_dim(0)->name(), "time") != 0) {
		_EXCEPTION1("Variable \"time\" does not have dimension \"time\" in file \"%s\"",
			strFilename.c_str());
	}

	// Calendar attribute
	NcAtt * attTimeCal = varTime->get_att("calendar");
	std::string strCalendar;
	if (attTimeCal == NULL) {
		if (fWarnOnMissingCalendar) {
			Announce("WARNING: Variable \"time\" is missing \"calendar\" attribute; assuming \"standard\"");
		}
		strCalendar = "standard";
	} else {
		strCalendar = attTimeCal->as_string(0);
	}
	Time::CalendarType eCalendarType =
		Time::CalendarTypeFromString(strCalendar);

	// Units attribute
	NcAtt * attTimeUnits = varTime->get_att("units");
	if (attTimeUnits == NULL) {
		_EXCEPTION1("Variable \"time\" is missing \"units\" attribute in file \"%s\"",
			strFilename.c_str());
	}
	std::string strTimeUnits = attTimeUnits->as_string(0);

	// Load in time data
	DataVector<int> vecTimeInt;
	DataVector<float> vecTimeFloat;
	DataVector<double> vecTimeDouble;
	DataVector<ncint64> vecTimeInt64;

	if (varTime->type() == ncInt) {
		vecTimeInt.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeInt[0]), dimTime->size());

	} else if (varTime->type() == ncFloat) {
		vecTimeFloat.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeFloat[0]), dimTime->size());

	} else if (varTime->type() == ncDouble) {
		vecTimeDouble.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeDouble[0]), dimTime->size());

	} else if (varTime->type() == ncInt64) {
		vecTimeInt64.Initialize(dimTime->size());
		varTime->set_cur((long)0);
		varTime->get(&(vecTimeInt64[0]), dimTime->size());

	} else {
		_EXCEPTION1("Variable \"time\" has invalid type "
			"(expected \"int\", \"int64\", \"float\" or \"double\")"
			" in file \"%s\"", strFilename.c_str());
	}

	for (int t = 0; t < dimTime->size(); t++) {
		Time time(eCalendarType);
		if (varTime->type() == ncInt) {
			time.FromCFCompliantUnitsOffsetInt(
				strTimeUnits,
				vecTimeInt[t]);

		} else if (varTime->type() == ncFloat) {
			time.FromCFCompliantUnitsOffsetDouble(
				strTimeUnits,
				static_cast<double>(vecTimeFloat[t]));

		} else if (varTime->type() == ncDouble) {
			time.FromCFCompliantUnitsOffsetDouble(
				strTimeUnits,
//.........这里部分代码省略.........

void CopyNcVar(
	NcFile & ncIn,
	NcFile & ncOut,
	const std::string & strVarName,
	bool fCopyAttributes,
	bool fCopyData
) {
	if (!ncIn.is_valid()) {
		_EXCEPTIONT("Invalid input file specified");
	}
	if (!ncOut.is_valid()) {
		_EXCEPTIONT("Invalid output file specified");
	}
	NcVar * var = ncIn.get_var(strVarName.c_str());
	if (var == NULL) {
		_EXCEPTION1("NetCDF file does not contain variable \"%s\"",
			strVarName.c_str());
	}

	NcVar * varOut;

	std::vector<NcDim *> dimOut;
	dimOut.resize(var->num_dims());

	std::vector<long> counts;
	counts.resize(var->num_dims());

	long nDataSize = 1;

	for (int d = 0; d < var->num_dims(); d++) {
		NcDim * dimA = var->get_dim(d);

		dimOut[d] = ncOut.get_dim(dimA->name());

		if (dimOut[d] == NULL) {
			if (dimA->is_unlimited()) {
				dimOut[d] = ncOut.add_dim(dimA->name());
			} else {
				dimOut[d] = ncOut.add_dim(dimA->name(), dimA->size());
			}

			if (dimOut[d] == NULL) {
				_EXCEPTION2("Failed to add dimension \"%s\" (%i) to file",
					dimA->name(), dimA->size());
			}
		}
		if (dimOut[d]->size() != dimA->size()) {
			if (dimA->is_unlimited() && !dimOut[d]->is_unlimited()) {
				_EXCEPTION2("Mismatch between input file dimension \"%s\" and "
					"output file dimension (UNLIMITED / %i)",
					dimA->name(), dimOut[d]->size());
			} else if (!dimA->is_unlimited() && dimOut[d]->is_unlimited()) {
				_EXCEPTION2("Mismatch between input file dimension \"%s\" and "
					"output file dimension (%i / UNLIMITED)",
					dimA->name(), dimA->size());
			} else if (!dimA->is_unlimited() && !dimOut[d]->is_unlimited()) {
				_EXCEPTION3("Mismatch between input file dimension \"%s\" and "
					"output file dimension (%i / %i)",
					dimA->name(), dimA->size(), dimOut[d]->size());
			}
		}

		counts[d] = dimA->size();
		nDataSize *= counts[d];
	}

	// ncByte / ncChar type
	if ((var->type() == ncByte) || (var->type() == ncChar)) {
		DataVector<char> data;
		data.Initialize(nDataSize);

		varOut =
			ncOut.add_var(
				var->name(), var->type(),
				dimOut.size(), (const NcDim**)&(dimOut[0]));

		if (varOut == NULL) {
			_EXCEPTION1("Cannot create variable \"%s\"", var->name());
		}

		var->get(&(data[0]), &(counts[0]));
		varOut->put(&(data[0]), &(counts[0]));
	}

	// ncShort type
	if (var->type() == ncShort) {
		DataVector<short> data;
		data.Initialize(nDataSize);

		varOut =
			ncOut.add_var(
				var->name(), var->type(),
				dimOut.size(), (const NcDim**)&(dimOut[0]));

		if (varOut == NULL) {
			_EXCEPTION1("Cannot create variable \"%s\"", var->name());
		}

		if (fCopyData) {
			var->get(&(data[0]), &(counts[0]));
//.........这里部分代码省略.........

int main(int argc, char ** argv) {

	MPI_Init(&argc, &argv);

try {

	// Number of latitudes
	int nLat;

	// Number of longitudes
	int nLon;

	// Zonal wavenumber
	int nK;

	// Meridional power
	int nLpow;

	// Output file
	std::string strOutputFile;

	// Parse the command line
	BeginCommandLine()
		CommandLineInt(nLat, "lat", 40);
		CommandLineInt(nLon, "lon", 80);
		CommandLineString(strOutputFile, "out", "topo.nc");

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	// Generate longitude and latitude arrays
	AnnounceBanner();
	AnnounceStartBlock("Generating longitude and latitude arrays");

	DataVector<double> dLon;
	dLon.Initialize(nLon);

	Parameters param;
	param.GenerateLatituteArray(nLat);

	std::vector<double> & dLat = param.vecNode;

	double dDeltaLon = 2.0 * M_PI / static_cast<double>(nLon);
	for (int i = 0; i < nLon; i++) {
		dLon[i] = (static_cast<double>(i) + 0.5) * dDeltaLon;
	}

	AnnounceEndBlock("Done");

	// Open NetCDF output file
	AnnounceStartBlock("Writing to file");

	NcFile ncdf_out(strOutputFile.c_str(), NcFile::Replace);

	// Output coordinates
	NcDim * dimLat = ncdf_out.add_dim("lat", nLat);
	NcDim * dimLon = ncdf_out.add_dim("lon", nLon);

	NcVar * varLat = ncdf_out.add_var("lat", ncDouble, dimLat);
	varLat->set_cur((long)0);
	varLat->put(&(param.vecNode[0]), nLat);

	NcVar * varLon = ncdf_out.add_var("lon", ncDouble, dimLon);
	varLon->set_cur((long)0);
	varLon->put(&(dLon[0]), nLon);

	// Generate topography
	DataMatrix<double> dTopo;
	dTopo.Initialize(nLat, nLon);

	double dK = static_cast<double>(nK);
	double dLpow = static_cast<double>(nLpow);

	double dA = 6.37122e6;
	double dX = 500.0;
	double dLatM = 0.0;
	double dLonM = M_PI / 4.0;
	double dD = 5000.0;
	double dH0 = 1.0;
	double dXiM = 4000.0;

	for (int j = 0; j < nLat; j++) {
	for (int i = 0; i < nLon; i++) {

		// Great circle distance
		double dR = dA / dX * acos(sin(dLatM) * sin(dLat[j])
			+ cos(dLatM) * cos(dLat[j]) * cos(dLon[i] - dLonM));

		double dCosXi = 1.0; //cos(M_PI * dR / dXiM);

		dTopo[j][i] = dH0 * exp(- dR * dR / (dD * dD))
			* dCosXi * dCosXi;
	}
	}

	// Write topography
	NcVar * varZs = ncdf_out.add_var("Zs", ncDouble, dimLat, dimLon);
	varZs->set_cur(0, 0);
	varZs->put(&(dTopo[0][0]), nLat, nLon);

//.........这里部分代码省略.........

int main(int argc, char ** argv) {

try {

	// Parameters
	Parameters param;

	// Output filename
	std::string strOutputFile;

	// Horizontal minimum wave number
	int nKmin;

	// Horizontal maximum wave number
	int nKmax;

	// Parse the command line
	BeginCommandLine()
		CommandLineInt(param.nPhiElements, "n", 40);
		CommandLineInt(nKmin, "kmin", 1);
		CommandLineInt(nKmax, "kmax", 20);
		CommandLineDouble(param.dXscale, "X", 1.0);
		CommandLineDouble(param.dT0, "T0", 300.0);
		CommandLineDouble(param.dU0, "U0", 20.0);
		CommandLineDouble(param.dG, "G", 9.80616);
		CommandLineDouble(param.dOmega, "omega", 7.29212e-5);
		CommandLineDouble(param.dGamma, "gamma", 1.4);
		CommandLineString(strOutputFile, "out", "wave.nc");

		ParseCommandLine(argc, argv);
	EndCommandLine(argv)

	AnnounceBanner();

	// Generate latitude values
	param.GenerateLatituteArray(param.nPhiElements);

	// Open NetCDF file
	NcFile ncdf_out(strOutputFile.c_str(), NcFile::Replace);

	NcDim *dimK = ncdf_out.add_dim("k", nKmax - nKmin + 1);
	NcDim *dimLat = ncdf_out.add_dim("lat", param.nPhiElements);
	NcDim *dimEig = ncdf_out.add_dim("eig", param.nPhiElements);

	// Write parameters and latitudes to file
	param.WriteToNcFile(ncdf_out, dimLat, dimLatS);

	// Wave numbers 
	NcVar *varK = ncdf_out.add_var("k", ncInt, dimK);

	DataVector<int> vecK;
	vecK.Initialize(nKmax - nKmin + 1);
	for (int nK = nKmin; nK <= nKmax; nK++) { 
		vecK[nK - nKmin] = nK;
	}

	varK->set_cur((long)0);
	varK->put(vecK, nKmax - nKmin + 1);

	// Eigenvalues
	NcVar *varMR = ncdf_out.add_var("mR", ncDouble, dimK, dimEig);
	NcVar *varMI = ncdf_out.add_var("mI", ncDouble, dimK, dimEig);

	NcVar *varUR = ncdf_out.add_var("uR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varUI = ncdf_out.add_var("uI", ncDouble, dimK, dimEig, dimLat);

	NcVar *varVR = ncdf_out.add_var("vR", ncDouble, dimK, dimEig, dimLatS);
	NcVar *varVI = ncdf_out.add_var("vI", ncDouble, dimK, dimEig, dimLatS);

	NcVar *varPR = ncdf_out.add_var("pR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varPI = ncdf_out.add_var("pI", ncDouble, dimK, dimEig, dimLat);

	NcVar *varWR = ncdf_out.add_var("wR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varWI = ncdf_out.add_var("wI", ncDouble, dimK, dimEig, dimLat);

	NcVar *varRhoR = ncdf_out.add_var("rhoR", ncDouble, dimK, dimEig, dimLat);
	NcVar *varRhoI = ncdf_out.add_var("rhoI", ncDouble, dimK, dimEig, dimLat);

	// Allocate temporary arrays
	DataVector<double> dUR;
	dUR.Initialize(param.nPhiElements);

	DataVector<double> dUI;
	dUI.Initialize(param.nPhiElements);

	DataVector<double> dVR;
	dVR.Initialize(param.nPhiElements-1);

	DataVector<double> dVI;
	dVI.Initialize(param.nPhiElements-1);

	DataVector<double> dPR;
	dPR.Initialize(param.nPhiElements);

	DataVector<double> dPI;
	dPI.Initialize(param.nPhiElements);

	DataVector<double> dWR;
	dWR.Initialize(param.nPhiElements);

//.........这里部分代码省略.........

void GridPatchCartesianGLL::InterpolateData(
	const DataVector<double> & dAlpha,
	const DataVector<double> & dBeta,
	const DataVector<int> & iPatch,
	DataType eDataType,
	DataLocation eDataLocation,
	bool fInterpAllVariables,
	DataMatrix3D<double> & dInterpData,
	bool fIncludeReferenceState,
	bool fConvertToPrimitive
) {
	if (dAlpha.GetRows() != dBeta.GetRows()) {
		_EXCEPTIONT("Point vectors must have equivalent length.");
	}

	// Vector for storage interpolated points
	DataVector<double> dAInterpCoeffs;
	dAInterpCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dBInterpCoeffs;
	dBInterpCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dADiffCoeffs;
	dADiffCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dBDiffCoeffs;
	dBDiffCoeffs.Initialize(m_nHorizontalOrder);

	DataVector<double> dAInterpPt;
	dAInterpPt.Initialize(m_nHorizontalOrder);

	// Element-wise grid spacing
	double dDeltaA = GetElementDeltaA();
	double dDeltaB = GetElementDeltaB();

	// Physical constants
	const PhysicalConstants & phys = m_grid.GetModel().GetPhysicalConstants();

	// Loop throught all points
	for (int i = 0; i < dAlpha.GetRows(); i++) {

		// Element index
		if (iPatch[i] != GetPatchIndex()) {
			continue;
		}

		// Verify point lies within domain of patch
		const double Eps = 1.0e-10;
		if ((dAlpha[i] < m_box.GetAEdge(m_box.GetAInteriorBegin()) - Eps) ||
			(dAlpha[i] > m_box.GetAEdge(m_box.GetAInteriorEnd()) + Eps) ||
			(dBeta[i] < m_box.GetBEdge(m_box.GetBInteriorBegin()) - Eps) ||
			(dBeta[i] > m_box.GetBEdge(m_box.GetBInteriorEnd()) + Eps)
		) {
			_EXCEPTIONT("Point out of range");
		}

		// Determine finite element index
		int iA =
			(dAlpha[i] - m_box.GetAEdge(m_box.GetAInteriorBegin())) / dDeltaA;

		int iB =
			(dBeta[i] - m_box.GetBEdge(m_box.GetBInteriorBegin())) / dDeltaB;

		// Bound the index within the element
		if (iA < 0) {
			iA = 0;
		}
		if (iA >= (m_box.GetAInteriorWidth() / m_nHorizontalOrder)) {
			iA = m_box.GetAInteriorWidth() / m_nHorizontalOrder - 1;
		}
		if (iB < 0) {
			iB = 0;
		}
		if (iB >= (m_box.GetBInteriorWidth() / m_nHorizontalOrder)) {
			iB = m_box.GetBInteriorWidth() / m_nHorizontalOrder - 1;
		}

		iA = m_box.GetHaloElements() + iA * m_nHorizontalOrder;
		iB = m_box.GetHaloElements() + iB * m_nHorizontalOrder;

		// Compute interpolation coefficients
		PolynomialInterp::LagrangianPolynomialCoeffs(
			m_nHorizontalOrder,
			&(m_box.GetAEdges()[iA]),
			dAInterpCoeffs,
			dAlpha[i]);

		PolynomialInterp::LagrangianPolynomialCoeffs(
			m_nHorizontalOrder,
			&(m_box.GetBEdges()[iB]),
			dBInterpCoeffs,
			dBeta[i]);

		int nComponents;
		int nRElements = m_grid.GetRElements();

		double ** pData2D;

		// State Data: Perform interpolation on all variables
		if (eDataType == DataType_State) {
//.........这里部分代码省略.........

void GridPatchCartesianGLL::EvaluateTestCase(
	const TestCase & test,
	const Time & time,
	int iDataIndex
) {
	// Initialize the data at each node
	if (m_datavecStateNode.size() == 0) {
		_EXCEPTIONT("InitializeData must be called before InitialConditions");
	}
	if (iDataIndex >= m_datavecStateNode.size()) {
		_EXCEPTIONT("Invalid iDataIndex (out of range)");
	}

	// Check dimensionality
	if ((m_grid.GetModel().GetEquationSet().GetDimensionality() == 2) &&
		(m_nVerticalOrder != 1)
	) {
		_EXCEPTIONT("VerticalOrder / Dimensionality mismatch:\n"
			"For 2D problems vertical order must be 1.");
	}

	// Evaluate topography
	EvaluateTopography(test);

	// Physical constants
	const PhysicalConstants & phys = m_grid.GetModel().GetPhysicalConstants();
	
	// Initialize the topography at each node
	for (int i = 0; i < m_box.GetATotalWidth(); i++) {
	for (int j = 0; j < m_box.GetBTotalWidth(); j++) {
		m_dataTopography[i][j] =
			test.EvaluateTopography(
				phys,
				m_dataLon[i][j],
				m_dataLat[i][j]);

		if (m_dataTopography[i][j] >= m_grid.GetZtop()) {
			_EXCEPTIONT("TestCase topography exceeds model top.");
		}

		// Gal-Chen and Sommerville vertical coordinate
		for (int k = 0; k < m_grid.GetRElements(); k++) {
			m_dataZLevels[k][i][j] =
				m_dataTopography[i][j]
					+ m_grid.GetREtaLevel(k)
						* (m_grid.GetZtop() - m_dataTopography[i][j]);
		}
		for (int k = 0; k <= m_grid.GetRElements(); k++) {
			m_dataZInterfaces[k][i][j] =
				m_dataTopography[i][j]
					+ m_grid.GetREtaInterface(k)
						* (m_grid.GetZtop() - m_dataTopography[i][j]);
		}

/*
		// Schar Exponential Decay vertical coordinate
		for (int k = 0; k < m_grid.GetRElements(); k++) {
			m_dataZLevels[k][i][j] = m_grid.GetZtop() * m_grid.GetREtaLevel(k) + 
			m_dataTopography[i][j] * sinh(m_grid.GetZtop() * (1.0 - m_grid.GetREtaLevel(k)) / m_dSL) / 
			sinh(m_grid.GetZtop() / m_dSL);
		}
		for (int k = 0; k <= m_grid.GetRElements(); k++) {
			m_dataZInterfaces[k][i][j] = m_grid.GetZtop() * m_grid.GetREtaInterface(k) + 
			m_dataTopography[i][j] * sinh(m_grid.GetZtop() * (1.0 - m_grid.GetREtaInterface(k)) / m_dSL) / 
			sinh(m_grid.GetZtop() / m_dSL);
		}
*/
	}
	}

	// Initialize the Rayleigh friction strength at each node
	if (test.HasRayleighFriction()) {
		for (int i = 0; i < m_box.GetATotalWidth(); i++) {
		for (int j = 0; j < m_box.GetBTotalWidth(); j++) {
			for (int k = 0; k < m_grid.GetRElements(); k++) {
				m_dataRayleighStrengthNode[k][i][j] =
					test.EvaluateRayleighStrength(
						m_dataZLevels[k][i][j],
						m_dataLon[i][j],
						m_dataLat[i][j]);
			}
			for (int k = 0; k < m_grid.GetRElements(); k++) {
				m_dataRayleighStrengthREdge[k][i][j] =
					test.EvaluateRayleighStrength(
						m_dataZInterfaces[k][i][j],
						m_dataLon[i][j],
						m_dataLat[i][j]);
			}
		}
		}
	}

	// Buffer vector for storing pointwise states
	const EquationSet & eqns = m_grid.GetModel().GetEquationSet();

	int nComponents = eqns.GetComponents();
	int nTracers = eqns.GetTracers();

	DataVector<double> dPointwiseState;
	dPointwiseState.Initialize(nComponents);
//.........这里部分代码省略.........