Commit 50d10e74 authored by Kirill Terekhov's avatar Kirill Terekhov

Added functionality and fixed issues.

Fixes for GCC.

Moved functionality to operate with dense matrices as INMOST service.

Added possibility to exchange data of type DATA_VARIABLE in parallel.

Added possibility to include external content in xml format.

Moved overloaded operators for variables from namespace INMOST into
global scope to prevent possible ambiguities in compilers.

Added perturbation algorithm for BCGS solver.
parent d03993fe
......@@ -240,6 +240,7 @@ set(INMOST_INSTALL_HEADERS Source/Headers/inmost.h
Source/Headers/inmost_autodiff.h
Source/Headers/inmost_common.h
Source/Headers/inmost_data.h
Source/Headers/inmost_dense.h
Source/Headers/inmost_expression.h
Source/Headers/inmost_mesh.h
Source/Headers/inmost_nonlinear.h
......@@ -268,6 +269,7 @@ set_property(TARGET inmost PROPERTY PUBLIC_HEADER
"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_autodiff.h"
"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_common.h"
"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_data.h"
"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_dense.h"
"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_expression.h"
"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_mesh.h"
"${PROJECT_SOURCE_DIR}/Source/Headers/inmost_nonlinear.h"
......
project(ADMFD)
set(SOURCE main.cpp matrix.hpp)
set(SOURCE main.cpp)
add_executable(ADMFD ${SOURCE})
......
......@@ -4,7 +4,6 @@
#include <string.h>
#include "inmost.h"
#include "matrix.hpp"
using namespace INMOST;
#ifndef M_PI
......@@ -184,8 +183,8 @@ int main(int argc,char ** argv)
rMatrix nKL(3,1); //normal to face
rMatrix lKs(3,1); //reminder of the co-normal vector Kn in cell K when projection to normal is substracted, i.e. (Kn - n n.Kn)
rMatrix lLs(3,1); //reminder of the co-normal vector Kn in cell L when projection to normal is substracted, i.e. (Kn - n n.Kn)
rMatrix KK = FromTensor(cK->RealArray(tag_K)).Transpose(); //diffusion tensor in cell K
rMatrix KL = FromTensor(cL->RealArray(tag_K)).Transpose(); //diffusion tensor in cell L
rMatrix KK = rMatrix::FromTensor(cK->RealArray(tag_K).data(),cK->RealArray(tag_K).size()).Transpose(); //diffusion tensor in cell K
rMatrix KL = rMatrix::FromTensor(cL->RealArray(tag_K).data(),cL->RealArray(tag_K).size()).Transpose(); //diffusion tensor in cell L
cK->Centroid(xK.data()); //retrive center of cell K
cL->Centroid(xL.data()); //retrive center of cell L
......@@ -237,7 +236,7 @@ int main(int argc,char ** argv)
real aF, vP = cell->Volume(); //area of the face
cell->Centroid(xP); //obtain cell center
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = faces.size(), k; //number of faces;
enumerator NF = (enumerator)faces.size(), k; //number of faces;
rMatrix IP(NF,NF), N(NF,3), R(NF,3); //matrices for inner product
......@@ -261,8 +260,8 @@ int main(int argc,char ** argv)
} //end of loop over faces
//inner product definition from Corollary 4.2, "Mimetic scalar products of differential forms" by Brezzi et al
IP = R*(R.Transpose()*N).Invert()*R.Transpose(); // Consistency part
IP += (rMatrix::Unit(NF) - N*(N.Transpose()*N).Invert()*N.Transpose())*(2.0/static_cast<real>(NF)*IP.Trace()); //Stability part
IP = R*(R.Transpose()*N).Invert().first*R.Transpose(); // Consistency part
IP += (rMatrix::Unit(NF) - N*(N.Transpose()*N).Invert().first*N.Transpose())*(2.0/static_cast<real>(NF)*IP.Trace()); //Stability part
//assert(IP.isSymmetric()); //test positive definitiness as well!
/*
......@@ -288,9 +287,9 @@ int main(int argc,char ** argv)
real nF[3]; //normal to the face
cell->Centroid(xP);
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = faces.size(), k; //number of faces;
enumerator NF = (enumerator)faces.size(), k; //number of faces;
dynarray<real,64> NG(NF,0), G(NF,0); // count how many times gradient was calculated for face
rMatrix K = FromTensor(cell->RealArrayDF(tag_K)); //get permeability for the cell
rMatrix K = rMatrix::FromTensor(cell->RealArrayDF(tag_K).data(),cell->RealArrayDF(tag_K).size()); //get permeability for the cell
rMatrix GRAD(NF,NF), NK(NF,3), R(NF,3); //big gradient matrix, co-normals, directions
rMatrix GRADF(NF,NF), NKF(NF,3), RF(NF,3); //corresponding matrices with masked negative values when multiplied by individual co-normals
rMatrix GRADFstab(NF,NF);
......@@ -305,18 +304,18 @@ int main(int argc,char ** argv)
R(k,1) = (yF[1]-xP[1]);
R(k,2) = (yF[2]-xP[2]);
// assemble matrix of co-normals
rMatrix nK = FromVector(nF).Transpose()*K;
rMatrix nK = rMatrix::FromVector(nF,3).Transpose()*K;
NK(k,0) = nK(0,0);
NK(k,1) = nK(0,1);
NK(k,2) = nK(0,2);
k++;
} //end of loop over faces
GRAD = NK*(NK.Transpose()*R).Invert()*NK.Transpose(); //stability part
GRAD = NK*(NK.Transpose()*R).Invert().first*NK.Transpose(); //stability part
//std::cout << "consistency" << std::endl;
//GRADF.Print();
GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert()*R.Transpose())*(2.0/NF*GRADF.Trace());
GRAD += (rMatrix::Unit(NF) - R*(R.Transpose()*R).Invert().first*R.Transpose())*(2.0/NF*GRADF.Trace());
/*
std::cout << "stability" << std::endl;
......@@ -400,7 +399,7 @@ int main(int argc,char ** argv)
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell)
{
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = faces.size(), k = 0;
enumerator NF = (enumerator)faces.size(), k = 0;
Cell cK = cell->self();
real gK = cK->Real(tag_F0), gL; //force on current cell and on adjacent cell
rMatrix M(cK->RealArrayDF(tag_M).data(),NF,NF); //inner product matrix
......@@ -475,7 +474,7 @@ int main(int argc,char ** argv)
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell) //loop over cells
{
ElementArray<Face> faces = cell->getFaces(); //obtain faces of the cell
enumerator NF = faces.size(), k = 0;
enumerator NF = (enumerator)faces.size(), k = 0;
Cell cK = cell->self();
variable pK = P(cK), pL; //pressure for back and front cells
rMatrix GRAD(cK->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
......@@ -509,7 +508,7 @@ int main(int argc,char ** argv)
for(Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell) //loop over cells
{
ElementArray<Face> facesK = cell->getFaces(), facesL; //obtain faces of the cell
enumerator NF = facesK.size(), k, l;
enumerator NF = (enumerator)facesK.size(), k, l;
Cell cK = cell->self();
variable pK = P(cK), pL; //pressure for back and front cells
rMatrix GRAD(cK->RealArrayDV(tag_W).data(),NF,NF); //Matrix for gradient
......
#ifndef __MATRIX_H
#define __MATRIX_H
#include "inmost.h"
using namespace INMOST;
// Matrix with n columns and m rows
// __m__
// | |
// n| |
// |_____|
//
// todo:
// 1. expression templates for operations
// 2. (ok) template matrix type for AD variables
template<typename Var>
class Matrix
{
public:
typedef unsigned enumerator;
protected:
array<Var> space;
enumerator n, m;
public:
static Matrix Unit(enumerator pn)
{
Matrix ret(pn,pn);
for(enumerator i = 0; i < pn; ++i)
{
for(enumerator j = 0; j < pn; ++j)
{
if( i == j )
ret(i,j) = 1;
else
ret(i,j) = 0;
}
}
return ret;
}
Matrix(Var * pspace, enumerator pn, enumerator pm) : space(pspace,pspace+pn*pm), n(pn), m(pm) {}
Matrix(enumerator pn, enumerator pm) : space(pn*pm), n(pn), m(pm) {}
Matrix(const Matrix & other) : space(other.n*other.m), n(other.n), m(other.m)
{
for(enumerator i = 0; i < n*m; ++i)
space[i] = other.space[i];
}
~Matrix() {remove();}
Matrix & operator =(Matrix const & other)
{
remove();
space.resize(other.n*other.m);
for(enumerator i = 0; i < other.n*other.m; ++i)
space[i] = other.space[i];
n = other.n;
m = other.m;
return *this;
}
void remove() {}
// i is in [0,n] - row index
// j is in [0,m] - column index
Var & operator()(enumerator i, enumerator j)
{
assert(i >= 0 && i < n);
assert(j >= 0 && j < m);
assert(i*m+j < n*m); //overflow check?
return space[i*m+j];
}
Var operator()(enumerator i, enumerator j) const
{
assert(i >= 0 && i < n);
assert(j >= 0 && j < m);
assert(i*m+j < n*m); //overflow check?
return space[i*m+j];
}
Matrix operator-(const Matrix & other) const
{
assert(n == other.n);
assert(m == other.m);
Matrix ret(n,m); //check RVO
for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]-other.space[k];
return ret;
}
Matrix & operator-=(const Matrix & other)
{
assert(n == other.n);
assert(m == other.m);
for(enumerator k = 0; k < n*m; ++k) space[k] -= other.space[k];
return *this;
}
Matrix operator+(const Matrix & other) const
{
assert(n == other.n);
assert(m == other.m);
Matrix ret(n,m); //check RVO
for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]+other.space[k];
return ret;
}
Matrix & operator+=(const Matrix & other)
{
assert(n == other.n);
assert(m == other.m);
for(enumerator k = 0; k < n*m; ++k) space[k] += other.space[k];
return *this;
}
Matrix operator*(Var coef) const
{
Matrix ret(n,m); //check RVO
for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]*coef;
return ret;
}
Matrix & operator*=(Var coef)
{
for(enumerator k = 0; k < n*m; ++k) space[k] *= coef;
return *this;
}
Matrix operator/(Var coef) const
{
Matrix ret(n,m); //check RVO
for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]/coef;
return ret;
}
Matrix & operator/=(Var coef)
{
for(enumerator k = 0; k < n*m; ++k) space[k] /= coef;
return *this;
}
Matrix operator*(const Matrix & other) const
{
assert(m == other.n);
Matrix ret(n,other.m); //check RVO
for(enumerator i = 0; i < n; ++i) //loop rows
{
for(enumerator j = 0; j < other.m; ++j) //loop columns
{
Var tmp = 0.0;
for(enumerator k = 0; k < m; ++k)
tmp += (*this)(i,k)*other(k,j);
ret(i,j) = tmp;
}
}
return ret;
}
/*
template<typename InType>
Matrix<InType> operator*(const Matrix<InType> & other)
{
assert(m == other.n);
Matrix<InType> ret(n,other.m); //check RVO
for(enumerator i = 0; i < n; ++i) //loop rows
{
for(enumerator j = 0; j < other.m; ++j) //loop columns
{
InType tmp = 0.0;
for(enumerator k = 0; k < m; ++k)
tmp += (*this)(i,k)*other(k,j);
ret(i,j) = tmp;
}
}
return ret;
}
template<typename InType>
Matrix operator*(const Matrix<InType> & other)
{
assert(m == other.n);
Matrix ret(n,other.m); //check RVO
for(enumerator i = 0; i < n; ++i) //loop rows
{
for(enumerator j = 0; j < other.m; ++j) //loop columns
{
Var tmp = 0.0;
for(enumerator k = 0; k < m; ++k)
tmp += (*this)(i,k)*other(k,j);
ret(i,j) = tmp;
}
}
return ret;
}
*/
Matrix Transpose() const
{
Matrix ret(m,n);
for(enumerator i = 0; i < n; ++i)
{
for(enumerator j = 0; j < m; ++j)
{
ret(j,i) = (*this)(i,j);
}
}
return ret;
}
Matrix Invert() const
{
Matrix ret(m,n);
Matrix At = Transpose(); //m by n matrix
Matrix AtB = At*Unit(n); //m by n matrix
Matrix AtA = At*(*this); //m by m matrix
enumerator * order = new enumerator [m];
for(enumerator i = 0; i < m; ++i) order[i] = i;
for(enumerator i = 0; i < m; i++)
{
enumerator maxk = i, maxq = i, temp2;
Var max, temp;
max = fabs(AtA(maxk,maxq));
//Find best pivot
for(enumerator k = i; k < m; k++) // over rows
{
for(enumerator q = i; q < m; q++) // over columns
{
if( fabs(AtA(k,q)) > max )
{
max = fabs(AtA(k,q));
maxk = k;
maxq = q;
}
}
}
//Exchange rows
if( maxk != i )
{
for(enumerator q = 0; q < m; q++) // over columns of A
{
temp = AtA(maxk,q);
AtA(maxk,q) = AtA(i,q);
AtA(i,q) = temp;
}
//exchange rhs
for(enumerator q = 0; q < n; q++) // over columns of B
{
temp = AtB(maxk,q);
AtB(maxk,q) = AtB(i,q);
AtB(i,q) = temp;
}
}
//Exchange columns
if( maxq != i )
{
for(enumerator k = 0; k < m; k++) //over rows
{
temp = AtA(k,maxq);
AtA(k,maxq) = AtA(k,i);
AtA(k,i) = temp;
}
//remember order in sol
{
temp2 = order[maxq];
order[maxq] = order[i];
order[i] = temp2;
}
}
if( fabs(AtA(i,i)) < 1.0e-54 )
{
bool ok = true;
for(enumerator k = 0; k < n; k++) // over columns of B
{
if( fabs(AtB(i,k)/1.0e-54) > 1 )
{
ok = false;
break;
}
}
if( ok ) AtA(i,i) = AtA(i,i) < 0.0 ? - 1.0e-12 : 1.0e-12;
else assert(false && "Matrix is singular");
}
for(enumerator k = i+1; k < m; k++)
{
AtA(i,k) /= AtA(i,i);
AtA(k,i) /= AtA(i,i);
}
for(enumerator k = i+1; k < m; k++)
for(enumerator q = i+1; q < m; q++)
{
AtA(k,q) -= AtA(k,i) * AtA(i,i) * AtA(i,q);
}
for(enumerator k = 0; k < n; k++)
{
for(enumerator j = i+1; j < m; j++) //iterate over columns of L
{
AtB(j,k) -= AtB(i,k) * AtA(j,i);
}
AtB(i,k) /= AtA(i,i);
}
}
for(enumerator k = 0; k < n; k++)
{
for(enumerator i = m; i-- > 0; ) //iterate over rows of U
for(enumerator j = i+1; j < m; j++)
{
AtB(i,k) -= AtB(j,k) * AtA(i,j);
}
for(enumerator i = 0; i < m; i++)
ret(order[i],k) = AtB(i,k);
}
delete [] order;
return ret;
}
void Zero()
{
for(enumerator i = 0; i < n*m; ++i) space[i] = 0.0;
}
Var Trace() const
{
assert(n == m);
Var ret = 0.0;
for(enumerator i = 0; i < n; ++i) ret += (*this)(i,i);
return ret;
}
Var * data() {return space.data();}
const Var * data() const {return space.data();}
enumerator Rows() const {return n;}
enumerator Cols() const {return m;}
void Print(Storage::real threshold = 1.0e-10) const
{
for(enumerator k = 0; k < n; ++k)
{
for(enumerator l = 0; l < m; ++l)
{
if( fabs((*this)(k,l)) > threshold )
std::cout << (*this)(k,l);
else
std::cout << 0;
std::cout << " ";
}
std::cout << std::endl;
}
}
bool isSymmetric() const
{
if( n != m ) return false;
for(enumerator k = 0; k < n; ++k)
{
for(enumerator l = k+1; l < n; ++l)
if( fabs((*this)(k,l)-(*this)(l,k)) > 1.0e-7 )
return false;
}
return true;
}
Var DotProduct(const Matrix & other) const
{
assert(n == other.n);
assert(m == other.m);
Var ret = 0.0;
for(enumerator i = 0; i < n*m; ++i) ret += space[i]*other.space[i];
return ret;
}
};
typedef Matrix<Storage::real> rMatrix; //shortcut for real matrix
typedef Matrix<Storage::var> vMatrix; //shortcut for matrix with variations
// Multiplication of matrices with promotion
__INLINE vMatrix operator *(const rMatrix & A, const vMatrix & B)
{
assert(A.Cols() == B.Rows());
vMatrix ret(A.Rows(),B.Cols()); //check RVO
for(rMatrix::enumerator i = 0; i < A.Rows(); ++i) //loop rows
{
for(rMatrix::enumerator j = 0; j < B.Cols(); ++j) //loop columns
{
variable tmp = 0.0;
for(rMatrix::enumerator k = 0; k < A.Cols(); ++k)
tmp += A(i,k)*B(k,j);
ret(i,j) = tmp;
}
}
return ret;
}
// Multiplication of matrices with promotion
__INLINE vMatrix operator *(const vMatrix & A, const rMatrix & B)
{
assert(A.Cols() == B.Rows());
vMatrix ret(A.Rows(),B.Cols()); //check RVO
for(rMatrix::enumerator i = 0; i < A.Rows(); ++i) //loop rows
{
for(rMatrix::enumerator j = 0; j < B.Cols(); ++j) //loop columns
{
variable tmp = 0.0;
for(rMatrix::enumerator k = 0; k < A.Cols(); ++k)
tmp += A(i,k)*B(k,j);
ret(i,j) = tmp;
}
}
return ret;
}
//retrive permeability matrix from stored data
__INLINE rMatrix FromTensor(Storage::real_array & K)
{
rMatrix Kc(3,3);
switch(K.size())
{
case 1: //scalar permeability tensor
Kc.Zero();
Kc(0,0) = Kc(1,1) = Kc(2,2) = K[0];
case 3:
Kc.Zero(); //diagonal permeability tensor
Kc(0,0) = K[0]; //KXX
Kc(1,1) = K[1]; //KYY
Kc(2,2) = K[2]; //KZZ
break;
case 6: //symmetric permeability tensor
Kc(0,0) = K[0]; //KXX
Kc(0,1) = Kc(1,0) = K[1]; //KXY
Kc(0,2) = Kc(2,0) = K[2]; //KXZ
Kc(1,1) = K[3]; //KYY
Kc(1,2) = Kc(2,1) = K[4]; //KYZ
Kc(2,2) = K[5]; //KZZ
break;
case 9: //full permeability tensor
Kc(0,0) = K[0]; //KXX
Kc(0,1) = K[1]; //KXY
Kc(0,2) = K[2]; //KXZ
Kc(1,0) = K[3]; //KYX
Kc(1,1) = K[4]; //KYY
Kc(1,2) = K[5]; //KYZ
Kc(2,0) = K[6]; //KZX
Kc(2,1) = K[7]; //KZY
Kc(2,2) = K[8]; //KZZ
break;
}
return Kc;
}
//retrive vector in matrix form from array
__INLINE rMatrix FromVector(Storage::real n[3])
{
return rMatrix(n,3,1);
}
//create diagonal matrix from array
__INLINE rMatrix FromDiagonal(Storage::real * r, rMatrix::enumerator size)
{
rMatrix ret(size,size);
ret.Zero();
for(rMatrix::enumerator k = 0; k < size; ++k) ret(k,k) = r[k];
return ret;
}
//create diagonal matrix from array of inversed values
__INLINE rMatrix FromDiagonalInverse(Storage::real * r, rMatrix::enumerator size)
{
rMatrix ret(size,size);
ret.Zero();
for(rMatrix::enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k];
return ret;
}
#endif
\ No newline at end of file
......@@ -201,15 +201,16 @@ namespace INMOST
INMOST_DATA_REAL_TYPE Automatizator::Derivative(expr & var, const Storage & e, Sparse::Row & out, Storage::real multiply, void * user_data)
{
Sparse::RowMerger & m = GetMerger();
INMOST_DATA_REAL_TYPE ret;
var.current_stencil.resize(1);
var.current_stencil[0] = stencil_pair(e->GetHandle(),1.0);
var.resize_for_stencil();
ret = EvaluateSub(var,0,ENUMUNDEF,user_data);
merger.PushRow(1.0,out);
DerivativeFill(var, 0, ENUMUNDEF, merger, multiply, user_data);
merger.RetriveRow(out);
merger.Clear();
m.PushRow(1.0,out);
DerivativeFill(var, 0, ENUMUNDEF, m, multiply, user_data);
m.RetriveRow(out);
m.Clear();
return ret*multiply;
}
#else
......
......@@ -11,7 +11,9 @@ namespace INMOST
case DATA_INTEGER: return "INTEGER";
case DATA_BULK: return "BULK";
case DATA_REFERENCE: return "REFERENCE";
#if defined(USE_AUTODIFF)
case DATA_VARIABLE: return "VARIABLE";
#endif
}
return "UNKNOWN";
}
......@@ -25,7 +27,9 @@ namespace INMOST
case DATA_INTEGER: return sizeof(INMOST_DATA_INTEGER_TYPE);
case DATA_REAL: return sizeof(INMOST_DATA_REAL_TYPE);
case DATA_REFERENCE: return sizeof(HandleType);
#if defined(USE_AUTODIFF)
case DATA_VARIABLE: return sizeof(variable);
#endif
}
return 0;
}
......@@ -39,7 +43,9 @@ namespace INMOST
case DATA_INTEGER: return sizeof(inner_integer_array);
case DATA_BULK: return sizeof(inner_bulk_array);
case DATA_REFERENCE: return sizeof(inner_reference_array);
#if defined(USE_AUTODIFF)
case DATA_VARIABLE: return sizeof(inner_variable_array);
#endif
}
return 0;
}
......@@ -73,13 +79,17 @@ namespace INMOST
else if( type == DATA_INTEGER ) new (adata) inner_integer_array (*static_cast<const inner_integer_array * >(bdata));
else if( type == DATA_BULK ) new (adata) inner_bulk_array (*static_cast<const inner_bulk_array * >(bdata));
else if( type == DATA_REFERENCE ) new (adata) inner_reference_array(*static_cast<const inner_reference_array * >(bdata));
#if defined(USE_AUTODIFF)
else if( type == DATA_VARIABLE ) new (adata) inner_variable_array (*static_cast<const inner_variable_array * >(bdata));
#endif
}
#if defined(USE_AUTODIFF)
else if( type == DATA_VARIABLE ) //have to call constructor
{
for(INMOST_DATA_ENUM_TYPE k = 0; k < data_size; ++k)
new (static_cast<variable *>(adata)+k) variable(*(static_cast<const variable *>(bdata)+k));
}
#endif
else // fixed size array
memcpy(adata,bdata,data_size*bytes);
}
......@@ -109,11 +119,13 @@ namespace INMOST
(*static_cast<inner_reference_array *> (adata)).~inner_reference_array();
new (adata) inner_reference_array();
}
#if defined(USE_AUTODIFF)
else if( type == DATA_VARIABLE )
{
(*static_cast<inner_variable_array *> (adata)).~inner_variable_array();
new (adata) inner_variable_array();
}
#endif
}
}
......@@ -177,7 +189,12 @@ namespace INMOST
case DATA_REAL: mem->bulk_data_type = INMOST_MPI_DATA_REAL_TYPE; break;
case DATA_INTEGER: mem->bulk_data_type = INMOST_MPI_DATA_INTEGER_TYPE; break;
case DATA_REFERENCE: mem->bulk_data_type = INMOST_MPI_DATA_ENUM_TYPE; break;
case DATA_VARIABLE: mem->bulk_data_type = INMOST_MPI_DATA_REAL_TYPE; break; //mixture of types, should not get here
#if defined(USE_AUTODIFF)
case DATA_VARIABLE:
if( !Sparse::HaveRowEntryType() ) Sparse::CreateRowEntryType();
mem->bulk_data_type = Sparse::GetRowEntryType();
break;
#endif
}
mem->bytes_size = DataTypeBytesSize(mem->dtype);
if(mem->size == ENUMUNDEF )
......@@ -481,7 +498,9 @@ namespace INMOST
case DATA_INTEGER: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast<void *>(&arr[it]); if( p != NULL ) (*static_cast<inner_integer_array *>( p )).~inner_integer_array(); } break;
case DATA_BULK: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast<void *>(&arr[it]); if( p != NULL ) (*static_cast<inner_bulk_array *>( p )).~inner_bulk_array(); } break;
case DATA_REFERENCE: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast<void *>(&arr[it]); if( p != NULL ) (*static_cast<inner_reference_array *>( p )).~inner_reference_array();} break;
#if defined(USE_AUTODIFF)
case DATA_VARIABLE: for(INMOST_DATA_ENUM_TYPE it = new_size; it < old_size; ++it) {void * p = static_cast<void *>(&arr[it]); if( p != NULL ) (*static_cast<inner_variable_array *>( p )).~inner_variable_array(); } break;
#endif
}
}
}
......@@ -499,7 +518,9 @@ namespace INMOST
case DATA_INTEGER: for(INMOST_DATA_ENUM_TYPE it = old_size; it < new_size; ++it) new ( &arr[it] ) inner_integer_array(); break;
case DATA_BULK: for(INMOST_DATA_ENUM_TYPE it = old_size; it < new_size; ++it) new ( &arr[it] ) inner_bulk_array(); break;