Немного о тетраэдальных сетках. Разбиение границы. Минусы тетраэдальной сетки.

Фактически, этот пост - продолжение этого.
Остановились мы на том, что построили тетраэдальную сетку из кубической. Следующая стадия - построить разбиение границы. Главная проблема в том, что разбиение границы должно быть согласованно с разбиением самой области - то есть, если у нас есть треугольник на границе, то этот треугольник полностью принадлежит одному из тетраэдров. Естественно есть множество способов сделать это. Простейший - берём прямоугольник разбиваем его на 4 треугольника (по 3 вершины разными способами), далее пробегаем по всем тетраэдрам и проверяем каждый тетраэдр на наличие этого треугольника. Проблема очевидна - на больших сетках это будет занимать гигантское количество времени.
Теперь если подумать о более рациональных способах, то можно сделать так: хранить отдельно параллелепипеды, находящиеся у границы сетки. Тогда сначала можно будет сделать разбиение внутренних, а затем внешних вместе с границей - тогда пробегать по тетраэдрам не придётся вообще(ну даже если придётся, то для каждого всего по 6 тетраэдрам, а не по всем, что есть - уже существенный плюс). Ну либо наоборот, порядок не важен.

Приведённые мой алгоритмы удобны в плане построения сеток для тестовых задач, либо для задач где надо сравнить решения на кубической и тетраэдальной сетки, однако построение "чисто" тераэдальной сетки даёт бóльшие возможности. С другой стороны это создаёт проблемы. Одна из главной проблема - это нормали. В вычислении нормалей проблем нет - векторное произведение сторон треугольника даст нам эту нормаль, главная проблема в определении ориентации нормали - её направления. Если с параллелепипедами всё просто - нормали для каждого параллелепипеда одинаковые(например, совпадают с координатными ортами), то для тетраэдров они могут быть направлены как угодно. И, если требуется ориентация нормали(что обычно так и есть), то надо для каждого тетраэдра разворачивать все нормали как надо, например, по направлению из тетраэдра, либо как-то хитро нумеровать сетку чтобы нормали расположились как надо. И то и другое требуют каких-никаких затрат.
Ещё один минус связан с размером тетраэдров в том или ином направлении(а возможно во всех сразу). Заключается в том, что при автоматическом построении сетки могут получаться очень маленькие тетраэдры. В зависимости от целей построения сетки это может не играть роли, а может быть минусом. Для графики, например, при аппроксимации некоторой модели тетраэдрами для отрисовки если маленький тетраэдр получается где-то внутри области, то это может не играть роли вообще. Для вычислительных задач это более важный критерий и может влиять на качество всего решения в целом.

Вещи, написанные мной в этих двух постах, в какой-то мере могут быть применены и к треугольным сеткам, для них алгоритмы будут формулироваться даже проще. Например, при построении из прямоугольной сетки треугольной разбиение границы остаётся тем же - перестраивать его не надо, да и само разбиение делается проще - по диагонали прямоугольника.

Немного о тетраэдальный сетках. Алгоритмы построения.

С начала кратко, о чём будет пост: о построение тетраэдальных сеток (сеток, элементы которых - тетраэдры), получение этих сеток из сеток на треугольных призмах и параллелепипедальных сеток (далее буду называть их кубическими -  это короче и проще выговаривать и читать).
Теперь то, что привело меня к этому: для курсового проекта требовалось строить тетраэдальные сетки, причём большие (в итоге получилось 132651 узлов и 750000 тетраэдров), что руками делать не очень быстро. Есть различные методы построения, но я выбрал достаточно простой (форма области была не важна) - строится кубическая сетка, затем она разбивается на тетраэдры. И столкнулся с некоторыми проблемами - так, как я сначала  разбивал сетку - она получала не комфорной, ну и первом делом я полез в гугл и, что самое печально я не нашёл того, чего искал. Возможно я конечно плохо искал, но всё равно. Поэтому путь это будет хотя бы тут.
Ну чтож, а теперь начнём с краткого обзора алгоритмов построения тетраэдальных сеток(и соответсвено треугольных, если смотреть двухмерное разбиение).
1. Алгоритм исчерпывания
Выходные данные: граничный фронт (множество треугольников - разбиение границы области)
Принцип работы алгоритма достаточно простой и понятный, берём треугольник из фронта, откладываем по нормали к нему точку на определённую длину, если нет пересечений с тетраэдрами то добавляем тетраэдр в разбиение, если есть пытаемся модифицировать его(переносим точку) и добавляем в разбиение. Исходный треугольник удаляем из фронта, полученные(из тетраэдра) добавляем в него. Заканчиваем когда больше не можем строить тетраэдры.
Плюсы  и минусы: явным минусом является необходимость начального фронта - его тоже ещё надо построить, а так же нахождение оптимальной для задачи стратегии выбора шага по нормали и модификации  тетраэдров. Плюсам является то, что не важна геометрия области.

2. Простой фронтальный метод
Это не название а разновидность метода. Алгоритм исчерпывания тоже является фронтальным методом, однако у них есть значительное отличие.
Входные данные: граничный фронт и внутренние узлы
Принцип работы алгоритма ещё проще, чем в предыдущем: берём треугольник из фронта, добавляем к нему узел и  получаем тетраэдр. Новые треугольники идут в фронт.
Плюсы и минусы: опять-таки входные данные, но тут их ещё больше требуется - кроме фронта нам нужны ещё внутренние точки. Плюсом является то, что с помощью данного алгоритма можно построить сетку, с определёнными особенностями, например тетраэдры вытянутые по оси z.

3. С помощью вложенных сеток
Входные данные: тетраэдальная сетка(возможно - один тетраэдр)
Принцип работы алгоритма: к исходной сетки строится вложенная(возможно несколько раз)
Плюсы и минусы: тут тоже исходные данные, если в предыдущих методах нам нужно было только разбиение границы, то здесь всей области. Ещё большим минусом является то, что если мы ходим построить сетку из большого числа тетраэдров она получится очень мелкой.

4. Построение из кубической сетки или сетки на треугольных призмах
Выходные данные: кубическая сетка или сетка на треугольных призмах
Принцип работы алгоритма: всё просто и понятно разбиваем исходные элементы на тетраэдры так, чтобы получилась комфорная сетка.
Плюсы и минусы: требование таких сеток большим минусом не является, особенно для кубической сетки - её строить легко, поэтому здесь метод даже выигрывает. Главным минусом является минус исходных сеток, по сравнению и "чисто тетраэдальными" сетками - тетраэдрами можно хорошо аппроксимировать более геометрически сложные области, чем параллелипипедами или призмами.

Конечно это далеко не все алгоритмы. Плюс есть алгоритмы, позволяющие улучшить качество тетраэдальных сеток, но их я рассматривать не буду. Перейду сразу к тому: как правильно разбить треугольную призму на тетраэдры, а потом как разбить куб.
 Пусть есть вот такая призма(не все невидимые линии правильно показаны, точнее показа правильно только одна из трёх, но не суть).

 Из неё получаются следующие три тетраэдра: (2, 4, 5, 6), (1, 2, 3, 6) и (1, 2, 4, 6).
Теперь перейдём к параллелепипедам: по суть мы разбиваем один параллелепипед на две призмы, а эти призмы на тетраэдры. Главная проблема - это правильная локальная нумерация призм друг по отношению к другу. Возьмём эту призму и добавим к ней ещё два узла, тем самым получим парллелепипед:

На кривости рисунка из-за копипаста обращать внимание не будем. Если взять симметричную нумерацию(7->2, 8->5), как с делал изначально сетка получится не комфорная, что в некоторых(даже думаю в большинстве) случаев будет проблемой. Правильно же локальную нумерацию во второй(левой) призме ввести следующим образом(глоб. -> лок.): 3 -> 1, 1 -> 2, 7 -> 3, 6 -> 4, 4 -> 5, 8 -> 6. И тогда, используя предыдущее разбиение мы получим наши заветные 6 тетраэдров: (2, 4, 5, 6), (1, 2, 3, 6), (1, 2, 4, 6), (1, 6, 4, 8), (3, 1, 7 8) и (3, 1, 6, 8).

Саму сетку мы разбили, однако осталось ещё одна проблема - разбиение границ. Об этом расскажу в следующий раз, а то пост и так получился каким-то длинным.

Visual Studio: небольшие странности с dll

Без большого объявления подробностей, надо было сделать две вещи:

1. Использовать dll
2. Компилировать с помощью bat-файла

 Экспортировалась из dll функция с простым прототипом:

extern "C" __declspec(dllexport) unsigned int doit(string& dir)

Возможно конечно я не правильно что-то делал(не extern "C", например), но если я компилировал и dll и использующую его программу в IDE любой из 2008, 2010 и 2012 Visual Studio всё работало прекрасно. Если же компиляция производилась, с помощью cl, то приложение вылетало. Два дня я запускал приложение и компиляцию на различных компьютерах и в итоге выяснил интересную вещь: если компилировать в IDE всё работает, если компилировать из-под консоли на Windows 7 работает только в случае, если используется VS 2012, в случае 2008 и 2010 - приложение вылетает. На Windows XP с точностью наоборот - на 2008 и 2010 (а даже 2005) работает, на 2012 - вылетает. Как только компилировать не пробовал. В итоге отчаялся и решил поэкспериментировать с самой прогой и изменил портотип на:

extern "C" __declspec(dllexport) unsigned int doit(char* windir)

И всё заработало! И в IDE и из-под консоли. Я не знаю из-за чего возникли эти странности, но вот так вот оно получилось.
Будьте аккуратнее с тем, что экспортируете.

Случайный поиск или на удивление быстрый поиск

Постановка задачи простая: найти глобальный минимум функции. Метод решения - постой случайных поиск.
После моих испытаний случайной сортировки я немного не доверял этому методу - за сколько времени он может найти решение с требуемой точностью. Оказалось, что при правильной оптимизации алгоритма это может произойти довольно быстро.

Агоритм решения задачи
Для начало определимся с входными данными для алгоритма: это заданная точность,  вероятность достижения этой точности, и прямоугольник, в котором ищим минимум.
Вообще фигурой может быть не прямоугольник, но это простейшая фигура и для неё легко построить нужное распределение, поэтому используется он.

Собственно алгоритм простой:

1. Вычисляем требуемое количество точек, чтобы с указанной вероятностью P найти минимум с указанной точностью;
2. Генерируем с равномерным распределением нужное количество точек и находим в какой из целевая функция принимает минимальное значение.

Алгоритм очень простой, теперь поговорим про две реализации, первая - которая приходит первой на ум, вторая - немного оптимизированный вариант. На счёт генерации - стандартный генератор C/С++ примерно воспроизводит равномерное распределение, поэтому использовать можно его.

Первая реализация
Всё просто на очередной итерации генерируем новую точку - если значение в ней меньше значения в текущем приближении, выбираем её за новое приближение, иначе приближение остаётся старым.
Плюсы: быстрота и компактность реализации.
Минусы: для достижения приличной точности в большой области требуется выполнить очень много вычислений(более миллиарда) и поэтому эта реализация очень медленная.
Ну а теперь моя версия, как улучшить реализацию алгоритма, не затратив слишком много памяти, но выйграв по времени.

Вторая реализация
Поскольку точки могут генерироваться независимости друг от друга, а так же вычисляться значение функции, то первое что пришло на ум - распараллелить. Осталось это сделать так, чтобы прога не занимала слишком много памяти(3 миллиарда double - это всё-таки уже много, если я не ошибся 22.4 Гб, если размер типа 8 байт). Поэтому реализовывалось всё так:
У нас есть n потоков, каждый поток должен посчитать по m точек, за раз он считает N точек(тоже параллельно), выбирает из них минимум и сравнивает с текущим приближением, далее, если надо меняет текущее приближение и генерирует ещё N точек и так пока не сгенериует их m штук. В конце у нас остаётся n чисел - минимум из каждого потока, их них выбирается уже глобальный минимум.
Число N можно высчитывать как-нибудь из сходя из общего числа точек на поток, но я просто бал фиксированно N = 10000 (были и другие варианты). Потому что меньших N (при N = 100) затраты на создание потоков были больше чем выгода.
Так что если нам надо хранить три числа(координаты точки и значение функции), то уйдёт всего 3*(n + n*N) ячеек памяти, что намного меньше чем 3*n*m (учитываем, что число потоков n на персональном компьютере невелико - от 2 до 16, примерно, а число m может быть несколько сотен миллионов, то выигрыш в памяти, по сравнению в всеобщим распараллеливанием виден).
Теперь зачем нужно делать это порциями по N штук, а не посчитать по m на потоке - по тем же самым причинам время и память, в рассматриваемом случае минимум из N ищется очень быстро(если ещё упорядочивать при добавлении то ещё быстрее) и памяти надо дополнительной, под 3*N, а не 3*m.
Плюсы: скорость - вот тут и было моё удивление, меньше за минуту на 4-х потоках обрабатывает до миллиарда точек.
Минусы: в реализации конечно не сильно труден, но значительно труднее чем первый метод.

Реальные замеры времени мне было лень делать, однако даже на глаз видно, что вторая реализация в разы быстрее.

Оценка времени работы алгоритмов. Виды роста затрат.

В последнее время было лень что-то писать. Очень лень. Но я собрался и решил написать о довольно лёгкой и в тоже время нужной теме - временные оценки для алгоритмов. Их цель - примерно оценить алгоритм по времени работы, в зависимости от входных данных и дать возможность сравнить алгоритмы, не имея их реализации. В прицепе, у всех алгоритмов есть особенности и если они[алгоритмы] легко и быстро реализуемы, то в конкретном случае сравнение на "практике" будет намного информативнее.
Об обозначениях: обычно используется O-натация, то есть обозначается как О(f(x)) или o(f(x)). Первое можно интерпретировать, как c*f(x), c>1, во втором случае 0<c<=1. Сама функция f(x) может быть, как непрерывной. так и дискретной, под x здесь понимается размер входных данных. Например x=n - длина массива, то функция будет дискретной, x=l - длина отрезка, функция будет непрерывной.
И теперь к самой теме: характеры роста, то есть вид функции f(x). Перечислять их буду в порядке возрастания временных затрат: сначала самые быстрые, затем все более медленные.

1. Полиномиальный рост
• f(x) = log(x) - логарифмический рост, считается одной из лучших оценок. Поскольку log(x) ~ x^a, 0<a<1 при больших x, то этот рост можно назвать минимальным из полиномиальных. Поскольку оценка определяется с точностью умножения на константу, то основание логарифма не важно.
• f(x) = x^a, a >1 - общий вид полиномиального роста, является приемлемой оценкой. Чем меньше показатель степени a, тем оценка лучше - времени меньше. Как было указанно выше  рост вида x*log(x) можно привести к этому типу роста
2. Экспоненциальный рост
Здесь небольшая путаница в терминологии,  при оценки алгоритмов часто любой не полиномиальный рост называют экспоненциальным.

• f(x) = e^x - экспоненциальный рост, рост в геометрической прогрессии. Является не очень хорошей оценкой для работы, однако обычно не критической. Если есть полиномиальный алгоритм, решающий эту задачу - лучше выбирать его.
• f(n) = n!, f(x) = Г(x) - факториальный рост, его спокойно можно включать в следующую группу, однако такие оценки тоже часто встречаются. Поясняю, почему здесь две функции - факториал достаточно известная оценка, но она определенна только для дискретных величин, для непрерывных общение факториала - гамма-функция.
3. Алгоритмы, от которых лучше бежать подальше
Врятле такие оценки вы когда-нибудь встретите, они здесь больше для сравнения - может быть и хуже.

• f(n) = sf(n) - суперфакториальный рост. sf(n) = 1! * 2! * .. * n! - произведение факториалов.  Растёт очень быстро.
• f(n) = (n!)^(n!) - я не думаю, что у такой зависимости есть название вообще. При увеличении n на единицу порядок порядка n увеличивается на единицу. Убунтовский калькулятор (9!)^(9!) на нетбуке уже считает оооочень долго(за 10 минут не посчитал), а виндосовский уже (7!)^(7!) не считает - переполнение.

И теперь, наглядное изображение написанного. Графики указанных функций.
Вот первые функции. 

Добавил (n!)^(n!) на график, только при n=1,2,3,4. Для сравнения: с (4!)^(4!) смогло сравняться только 32!. Ну а что дальше - понятно.

Так что, выбирая алгоритм будьте осторожны - одно и туже задачу, бывает, можно решить очень разными способами.

Стандартный random C++

Что-то вдруг захотелось посмотреть, как работает функция rand() в C++. Не алгоритм генерации, а распределение результатов. Ну и забавы ради посмотрел что получилось.
Результаты получились довольной скучные - распределение примерно равномерное. Как проводились тесты - создаётся массив, зануляется, затем 10000000 прибавляется единица к случайному элементу массива, ну и затем выводится относительная частота. Объяснение немного странное, поэтому в конце поста есть код программы.

1. Без перемешивания, 32768 элементов.
 
 Под перемешиванием я имею ввиду функцию srand





 Первые 100 элементов из этой последовательности:














Перемешивание в начале генерации послеовательности, 32768 элементов
 
 Аналогично первые 100:











Перемешивание перед генерацией каждого числа, 256 элементов
 Вообще на самом деле не лучшая идея, т.к. по сути каждый раз генератор сбрасывается. Зато получается хоть какая-то оригинальная картинка.

Код программы(для последнего случая):

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <sstream>
#include <time.h>

using namespace std;

string convertInt(int number){
   stringstream ss;
   ss << number;
   return ss.str();
}

int main(){
 string f_n;
 FILE* outp;
 unsigned long int *mas = new unsigned long int [32769];
 int rand_N = 10000000;
 
 for(int n = 2; n < 32769; n *= 2){

  for(int i = 0; i < n; i++)
   mas[i] = 0;

  f_n = convertInt(n);
  f_n = f_n + "_tsr.txt";

  outp = fopen(f_n.c_str(), "w");
  for(int i = 0; i < rand_N; i++){
   mas[rand()%n]++;
   srand(time(0));
  }

  for(int i = 0; i < n; i++)
   fprintf(outp, "%d\t%f\n", i, 1.0*mas[i]/rand_N);
  fclose(outp);
 }
 return 0;
}

Для тестирования использовалась Visual Studio 2008 с её компилятором. Не думаю, что на gcc отличия будут большие(в среднем конечно) и получится примерно та же картинка.