Добрый день, коллеги.
Предо мной встала задача: по заданным траекториям скважин в формате welltrack и файлом с KW ZCORN и COORD определить ячейки, через которые проходят скважины. Чтобы это сделеать - нужно сначала считать сабж. С COORD вроде бы разобрался худо бедно. А с ZCORN не совсем: как последовательность координт TVDSS считать в нужные контейнеры типов точек (под типами точек понимаю здесь различные 8 типов точек (образованных ZCORNами и COORDами) для некоторго рассматриваемого блока - например верхняя левая угловая, верхняя правая угловая и т.д.). Получается на входе у меня есть только DIMENS, и порядковый номер считываемой позиции ZCORN(какое по счету число типа float я уже считал). Из этих данных мне нужно однозначно понять: к какому блоку принадлежит эта точка и какой у неё тип.
Кто нибудь решал такую/подобную задачу? НА форуме порылся, но не нашел похожих ответов. Скорее всего EmptyEye13 занимался этим, полагаю.
Пишу пока что в Python.
Спасибо, друзья за помощь
Посмотрите File Formats Reference Manual - спецификация файлов открытая, а вот что такое формат welltrack - не знаю.
Формат welltrack это всего-навсего последовательность значений:
-----------------------------------------------
-- Depth Ref. : Vertical Datum (e.g. MSL)
-- Position Ref. : Cartographic (e.g. UTM)
-- Depth Unit : [m]
-- Position Unit : [m]
-- Type : Actual Trajectory
-- RKB Elevation : 173.220000
welltrack 'wellname'
X1 Y1 TVDSS1 MD1
X2 Y2 TVDSS2 MD2
По поводу File Formats Reference Manual - не напомните, где бы я его мог найти? Обязательно посмотрю
После инсталляции есть ссылка на титульной странице Bookshelf на все мануалы
public Cell GetCell(int X, int Y, int Z)
{
// Формат именования вершин в кубе.
// На первом месте либо T (top, верхняя грань), либо B (bottom, нижняя грань)
// далее N (north, северная, условный верх) либо S (south, южная, условный низ) грань
// и завершается W( west, западная, условное лево) либо E (east, восточное, условное право).
//Таким образом, трехбуквенным кодом обозначаются восемь вершин одной ячейки.
// Это распространенный подход.
Cell CELL = new Cell();
// Отметки глубин
CELL.TNW.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + 2 * X + 0)];
CELL.TNE.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + 2 * X + 1)];
CELL.TSW.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + 2 * X + Project.NX * 2 + 0)];
CELL.TSE.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + 2 * X + Project.NX * 2 + 1)];
CELL.BNW.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + Project.NX * Project.NY * 4 + 2 * X + 0)];
CELL.BNE.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + Project.NX * Project.NY * 4 + 2 * X + 1)];
CELL.BSW.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + Project.NX * Project.NY * 4 + 2 * X + Project.NX * 2 + 0)];
CELL.BSE.Z = Project.ZCORN[(ulong)(Z * Project.NX * Project.NY * 8 + Y * Project.NX * 4 + Project.NX * Project.NY * 4 + 2 * X + Project.NX * 2 + 1)];
// Направляющая линия от TNW до BNW
Vector3 TOP, BTM;
TOP.X = Project.COORD[(X + (Project.NX + 1) * Y) * 6 + 0];
TOP.Y = Project.COORD[(X + (Project.NX + 1) * Y) * 6 + 1];
TOP.Z = Project.COORD[(X + (Project.NX + 1) * Y) * 6 + 2];
BTM.X = Project.COORD[(X + (Project.NX + 1) * Y) * 6 + 3 + 0];
BTM.Y = Project.COORD[(X + (Project.NX + 1) * Y) * 6 + 3 + 1];
BTM.Z = Project.COORD[(X + (Project.NX + 1) * Y) * 6 + 3 + 2];
float FRAC = 0;
if (BTM.Z == TOP.Z) // нет наклона направляющей линии, значит координаты равны
{
CELL.TNW.X = TOP.X;
CELL.TNW.Y = TOP.Y;
CELL.BNW.X = BTM.X;
CELL.BNW.Y = BTM.Y;
}
else
{
FRAC = (CELL.TNW.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.TNW.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.TNW.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
FRAC = (CELL.BNW.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.BNW.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.BNW.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
}
// Направляющая линия от TNE до BNE
TOP.X = Project.COORD[((X + 1) + (Project.NX + 1) * Y) * 6 + 0];
TOP.Y = Project.COORD[((X + 1) + (Project.NX + 1) * Y) * 6 + 1];
TOP.Z = Project.COORD[((X + 1) + (Project.NX + 1) * Y) * 6 + 2];
BTM.X = Project.COORD[((X + 1) + (Project.NX + 1) * Y) * 6 + 3 + 0];
BTM.Y = Project.COORD[((X + 1) + (Project.NX + 1) * Y) * 6 + 3 + 1];
BTM.Z = Project.COORD[((X + 1) + (Project.NX + 1) * Y) * 6 + 3 + 2];
if (BTM.Z == TOP.Z) // нет наклона направляющей линии, значит координаты равны
{
CELL.TNE.X = TOP.X;
CELL.TNE.Y = TOP.Y;
CELL.BNE.X = BTM.X;
CELL.BNE.Y = BTM.Y;
}
else
{
FRAC = (CELL.TNE.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.TNE.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.TNE.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
FRAC = (CELL.BNE.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.BNE.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.BNE.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
}
// Направляющая линия от TSE до BSE
TOP.X = Project.COORD[((X + 1) + (Project.NX + 1) * (Y + 1)) * 6 + 0];
TOP.Y = Project.COORD[((X + 1) + (Project.NX + 1) * (Y + 1)) * 6 + 1];
TOP.Z = Project.COORD[((X + 1) + (Project.NX + 1) * (Y + 1)) * 6 + 2];
BTM.X = Project.COORD[((X + 1) + (Project.NX + 1) * (Y + 1)) * 6 + 3 + 0];
BTM.Y = Project.COORD[((X + 1) + (Project.NX + 1) * (Y + 1)) * 6 + 3 + 1];
BTM.Z = Project.COORD[((X + 1) + (Project.NX + 1) * (Y + 1)) * 6 + 3 + 2];
if (BTM.Z == TOP.Z) // нет наклона направляющей линии, значит координаты равны
{
CELL.TSE.X = TOP.X;
CELL.TSE.Y = TOP.Y;
CELL.BSE.X = BTM.X;
CELL.BSE.Y = BTM.Y;
}
else
{
FRAC = (CELL.TSE.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.TSE.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.TSE.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
FRAC = (CELL.BSE.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.BSE.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.BSE.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
}
// Направляющая линия от TSW до BSW
TOP.X = Project.COORD[(X + (Project.NX + 1) * (Y + 1)) * 6 + 0];
TOP.Y = Project.COORD[(X + (Project.NX + 1) * (Y + 1)) * 6 + 1];
TOP.Z = Project.COORD[(X + (Project.NX + 1) * (Y + 1)) * 6 + 2];
BTM.X = Project.COORD[(X + (Project.NX + 1) * (Y + 1)) * 6 + 3 + 0];
BTM.Y = Project.COORD[(X + (Project.NX + 1) * (Y + 1)) * 6 + 3 + 1];
BTM.Z = Project.COORD[(X + (Project.NX + 1) * (Y + 1)) * 6 + 3 + 2];
if (BTM.Z == TOP.Z) // нет наклона направляющей линии, значит координаты равны
{
CELL.TSW.X = TOP.X;
CELL.TSW.Y = TOP.Y;
CELL.BSW.X = BTM.X;
CELL.BSW.Y = BTM.Y;
}
else
{
FRAC = (CELL.TSW.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.TSW.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.TSW.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
FRAC = (CELL.BSW.Z - TOP.Z) / (BTM.Z - TOP.Z);
CELL.BSW.X = TOP.X + FRAC * (BTM.X - TOP.X);
CELL.BSW.Y = TOP.Y + FRAC * (BTM.Y - TOP.Y);
}
return CELL;
}
public struct Cell
{
public Vector3 TNW;
public Vector3 TNE;
public Vector3 TSW;
public Vector3 TSE;
public Vector3 BNW;
public Vector3 BNE;
public Vector3 BSW;
public Vector3 BSE;
}
Обычно нужно не по индексу ZCORN понять к какой ячейке он принадлежит, а наоборот, по индексам ячейки найти высоты ее углов. Допустим, у нас для сетки с размерностью nx,ny,nz массив ZCORN прочитан в массив в памяти (индексация с нуля) и нам нужно получить высоты углов ячейки имеющей индексы i,j,k. Тогда углы этой ячейки (a, b, c, d - верхняя грань; e, f, g, h - нижняя) можно достать из массива ZCORN так:
nx2 = nx * 2;
ny2 = ny * 2;
nxy2 = nx2 * ny2;
i--;
j--;
k--;
a = zcorn[ k * nxy2 + j * nx2 + i * 2 ];
b = zcorn[ k * nxy2 + j * nx2 + i * 2 + 1 ];
c = zcorn[ k * nxy2 + (j+1) * nx2 + i * 2 ];
d = zcorn[ k * nxy2 + (j+1) * nx2 + i * 2 + 1 ];
e = zcorn[ (k+1) * nxy2 + j * nx2 + i * 2 ];
f = zcorn[ (k+1) * nxy2 + j * nx2 + i * 2 + 1 ];
g = zcorn[ (k+1) * nxy2 + (j+1) * nx2 + i * 2 ];
h = zcorn[ (k+1) * nxy2 + (j+1) * nx2 + i * 2 + 1 ];
Если же действительно нужно получать индексы ячейки по позиции в массиве ZCORN, то приведенные выше несколько строк кода можно без труда обратить.
Хотя IMHO делать вычисления при чтении больших объемов данных - путь к сложному в понимании и поддержке коду. Гораздо проще сначала все считать в память, а потом уже (обычно уже совсем в другом куске кода) разбираться что и от какой ячейки где лежит.
Вторая часть вопроса, когда-то волновала меня.
При определении пересечения траектории скважины с гранью ячейки, грань ячейки описывается билинейной поверхностью, в координатах u, v. Найти пересечение точки и билинейной поверхности это не сложная задача. Найти же пересечение элемента траектории - увы, была сложнее чем я мог решить в то время.
Можно поступить проще, грань ячейки апроксимируется четырьмя треугольниками, через добавления центральной точки на грань. Найти центр можно из билинейной интерполяции при u=0.5, v=0.5. Найти пересечение треугольника и линии это задача намного проще.
Хотя честное решение существует :)
Формат файла COORD/ZCORN и как это работает есть в мануалах Eclipse. Я когда-то использовал данное описание чтобы наборот создать дата файл из Excel с наклонными поверхностями, документации оказалось достаточно.
По поводу пересечения траектории и грани ячеек это обычная геометрическая задача, я бы сделал так:
- разбить траекторию скважины на линейные сегменты, скажем 50 метров
- провести первый фильтр и посмотреть какие ячейки из модели вообще рядом с каждым сегментом, сравнить расстояния между вершинами ячеек и двумя точками траектории данного сегмента (для скорости так как ячеек может быть очень много)
- так у нас на каждый сегмент будет небольшой набор ячеек которые возможно пересекает траектория.
- потом для каждой грани из этих потенциальных кандидатов сделать проверку на пересечение, это обычная геометрическая задача пересечения поверхности и прямой.
- найти эту точку пересечения если существует и проверить если она ограничена гранью ячейки.
- получим все грани которые пересекает кривая, далее уже тривиально сказать какие ячейки.
Добрый день,
Я прикрепил свой старый код на питоне который выводил 3д сетку с использованием либы http://docs.enthought.com/mayavi/mayavi/
это старый код, я не уверен что он сразу сработает на последней версии библиотеки, но он короткий и возможно поможет разобраться
Интересно/ хотелось бы увидеть рисунок с углами ячейки (a, b, c, d - верхняя грань; e, f, g, h - нижняя) с нанесением ортов \vec{i}, \vec{j} , \vec{k}, ведь от того, как мы пронумеруем b[ - будет зависеть то, какая формула применится.
СПасибо
Спасибо, EmptyEye13
Не совсем понимаю для чего нужно разбивать на 4 треугольника (я так понимаю речь идет только про верхнюю и нижнюю грани), если можно найти уравнения двух плоскостей, составлюящих данные поверхности и искать пересечение траектории уже с ними. а потом проверять пространственными неравенствами (если определитель будет больше нуля) - принадлежит ли точка пересечения объему ячейки.
Спасибо, RomanK.
Ну это просто. Потому-что грань ячейки не является плоскостью. Я выше написал, что это билинейная поверхность.
Линейная поверхность описывается тремя точками, билинейная - четырьмя. Для упрощения, можно приблизительно описать билинейную поверхность разбивая её на четыре линейных поверхности (четыре треугольника).
Я так понял, пока разбирался в своё время.
Спасибо. Теперь мне тоже ясно стало.
Я таких ячеек не встречал, теоретически они возможны, но мне что-то подсказывает большинство пакетов по геомоделировнию не делает таких ячеек или все таки верхняя крышка очень близка к плоскости на реальных моделях что это не бросается в глаза и тогда стоит ли заморачиваться. Еще мне кажется там были какие-то ограничения для corner-point geometry, что это не свободные направляющие.
Интересно как пакеты считают объемы таких ячеек, какое приближение используют.
Если считать что на картинке с заголовком "b." изображена ячейка с плоской нижней (серой) гранью, то на ней углы имеют подписи: a - 91, b - 210, c - 162, d - 95, и углы e,f,g,h - это те что соответственно под ними на нижней грани. И ещё: слои всего - углов, ячеек, значений в ячейках и т.п. задаются не снизу-вверх, а сверху-вниз (стрелочка "value" должна смотреть в обратную сторону). \vec{i}, \vec{j} , \vec{k} это стрелки "col", "row" и "value в обратном направлении" .
тоже так думаю
Ячейки могут попадаться абсолютно разные (особенно если ничем не ограничивать фантазию геолога при создании сетки модели :) ). Однажды даже встречалась настолько кривая ячейка (напоминающая по форме винт самолета) что ее центр тяжести лежал вне самой ячейки! Самодельный симулятор на ней крешился, но Еclipse нормально обрабатывал.
Подозреваю, что наиболее профессиональный открытый исходный код на эту тему можно найти тут:
https://github.com/OPM/ResInsight
В нем, кстати, уже решена проблема вычисления пресечения траектории скважины с ячейками сетки модели ;)
Насчет объемов - насколько я понимаю считается определитель смешанного произведения векторов
Полностью согласен. Спасибо за описание легенды рисунка, кстати=)
Тогда они считаются неправильно если есть такие сложные ячейки. Хотя можно посчитать объемы по всем 4 тройкам векторов и усреднить, но это будет уже приближение.
Объем таких ячеек, считается через формулу Остроградского-Гаусса и это делается в любом пакете геологического моделирования, также симуляторы считают объем ячейки при инициализации. Не смог пока найти в OPM реализованного алгоритма.
Приведу интересный финт, как геологические пакеты рассчитывают усеченные ячейки, например когда ячейка делится поверхностью ВНК и надо подсчитать объем над и под ВНК. В RMS объем ячейки заполняется равномерно точками и далее подсчитывается количество точек в каждой области, например,
Объем под ВНК = Общий объем * (N точек под ВНК / N точек)
RomanK: "Не смог пока найти в OPM реализованного алгоритма."
https://github.com/OPM/opm-grid/blob/master/opm/grid/cornerpoint_grid.c (void compute_geometry(struct UnstructuredGrid *g))
и
https://github.com/OPM/opm-grid/blob/master/opm/grid/cpgpreprocess/geometry.c (compute_cell_geometry и compute_cell_geometry_3d) ?