Commit 9bbc4789 authored by user's avatar user
Browse files

initial commit 

	new file:   Makefile
	new file:   coefs.cpp
	new file:   coefs.h
	new file:   core.cpp
	new file:   core.h
	new file:   cs.h
	new file:   defs.h
	new file:   lbm.cpp
	new file:   lbm.h
	new file:   main.cpp
	new file:   taylor.cpp
	new file:   taylor.h
parent 819f29e9

Makefile

0 → 100644
+22 −0
Original line number Diff line number Diff line 
defaults: main



CXXFLAGS=-O2 -I. -std=c++11

LIBS=-lm -lginac 



OBJS = lbm.o core.o main.o coefs.o taylor.o

DEPS = *.h



main: $(OBJS) $(DEPS)

	g++ $(CXXFLAGS) $(OBJS) -o main  $(LIBS) 



lbm.o: lbm.cpp defs.h lbm.h

	$(CXX) -MD -MP $(CXXFLAGS2) -c -o $@ $<



%.o: %.cpp $(DEPS)

	$(CXX) -MD -MP $(CXXFLAGS) -c -o $@ $<



clean:

	@rm -rf main

	@rm -rf *.o

	@rm -rf *.d

coefs.cpp

0 → 100644
+95 −0
Original line number Diff line number Diff line 
#include "coefs.h"



void COEFS::fill(ex &e)

{

	if (basis.size()<=0) printf("error: COEFS: no basis %d\n", basis.size());

	else  

	{

		SC sc;

		sc.init(basis.size());

		sc.coef=(ex)0;

		sc.base=(ex)1;

		fill_recursive(e, sc, 0);  // musts be expanded!!!!

	}

}



bool COEFS::fill_recursive(ex &e, SC &sc, int level)

{

	if (level>basis.size()) return false;

	if (level==basis.size()) // last entry

	{

		ex E=expand(e);

		if (E!=0) //  non-zero coefficient

		{

			sc.coef=E;

			v.push_back(sc);

			return true;

		} else

		return false;

	}

	for (int ii=0;ii<order;ii++)

	{

		// hotfix to enforce ordering

		int i = (ii==order-1) ? 0 : ii+1;

		if (sc.order()+i<order)

		{

			ex c = coeff(e, basis[level], i);

			if (c!=0)

			{

				sc.exp[level] = i;

				ex oldbase=sc.base;

				sc.base *= pow(basis[level],i);

				fill_recursive(c, sc, level+1);

				sc.base = oldbase;

				sc.exp[level] = 0;

			}

		}

	}

	return true;

}



ex COEFS::check(ex &e)

{

	// check whether all coefficients were correctly read

	ex expr = 0;

	for (int i=0;i<v.size();i++)

	{

		expr += v[i].coef*v[i].base;

	}

	ex diff=expand(e-expr);

	if (diff!=0)

	{

		cout << "check: nonzero diff = $"<<endl << diff <<"$"<< endl<< endl;

		cout << "source: $"<<endl << e<<endl <<"$"<< endl<< endl;

		cout << "target: $"<<endl << expr<<endl <<"$"<< endl<< endl;

	}

	return diff;

}



bool COEFS::compare(SC &sc1, SC &sc2)

{

	if (sc1.N != sc2.N) return false;

	for (int i=0;i<sc1.N;i++)

		if (sc1.exp[i]!=sc2.exp[i]) return false;

	return true;

}



void COEFS::add(COEFS &CF2)

{

	for (int j=0;j<CF2.v.size();j++)

	{

		// find match for 

		int match=-1;

		for (int i=0;i<v.size();i++)

			if (compare(v[i],CF2.v[j])) {  match = j; break; }

		if (match<0)

			v.push_back(CF2.v[j]);

		else

			v[match].coef += CF2.v[j].coef;

	}

}



void COEFS::copy_basis(COEFS &CF2)

{

	for (int i=0;i<CF2.basis.size();i++) basis.push_back(CF2.basis[i]);

}

coefs.h

0 → 100644
+38 −0
Original line number Diff line number Diff line 
#pragma once



#include "defs.h"



const int MAXSC = 1000; // limit



struct SC 

{

	int N;

	int exp[MAXSC]; // exponents

	ex coef;

	ex base;

	int order() { int sum=0; for (int i=0;i<N;i++) sum+=exp[i]; return sum; }

	void init(int iN)

	{ 

		N=iN;

		if (N>=MAXSC) printf("error: MAXSC %d exceeded %d\n",MAXSC, N);

		for (int i=0;i<N;i++) exp[i]=0; 

	}

};





struct COEFS 

{

	std::vector<SC> v;

	std::vector<ex> basis;

	std::vector<ex> basis_tex;

	int order=6; // default

	void fill(ex &e);

	bool fill_recursive(ex &e, SC &sc, int level);

	ex check(ex &e);

	// adds another COEFS to this

	void add(COEFS &CF2);

	// copy basis

	void copy_basis(COEFS &CF2);

	// compare exp[] for two SCs 

	bool compare(SC&sc1, SC&sc2);

};

core.cpp

0 → 100644
+1323 −0

File added.

Preview size limit exceeded, changes collapsed.

core.h

0 → 100644
+118 −0
Original line number Diff line number Diff line 
//  .----------------.  .----------------.  .----------------.  .----------------.  .-----------------.

// | .--------------. || .--------------. || .--------------. || .--------------. || .--------------. |

// | |   _____      | || |   ______     | || | ____    ____ | || |      __      | || | ____  _____  | |

// | |  |_   _|     | || |  |_   _ \    | || ||_   \  /   _|| || |     /  \     | || ||_   \|_   _| | |

// | |    | |       | || |    | |_) |   | || |  |   \/   |  | || |    / /\ \    | || |  |   \ | |   | |

// | |    | |   _   | || |    |  __'.   | || |  | |\  /| |  | || |   / ____ \   | || |  | |\ \| |   | |

// | |   _| |__/ |  | || |   _| |__) |  | || | _| |_\/_| |_ | || | _/ /    \ \_ | || | _| |_\   |_  | |

// | |  |________|  | || |  |_______/   | || ||_____||_____|| || ||____|  |____|| || ||_____|\____| | |

// | |              | || |              | || |              | || |              | || |              | |

// | '--------------' || '--------------' || '--------------' || '--------------' || '--------------' |

//  '----------------'  '----------------'  '----------------'  '----------------'  '----------------' 

//

// https://mmg-gitlab.fjfi.cvut.cz/gitlab/lbm/lbman

//

// Written by Radek Fučík (fucik@fjfi.cvut.cz) FNSPE CTU in Prague, Czech Republic

// with a great help of Robert Straka (straka@metal.agh.edu.pl) AGH Krakow, Poland

//

// This program is free software: you can redistribute it and/or modify

// it under the terms of the GNU General Public License as published by

// the Free Software Foundation, either version 3 of the License, or

// (at your option) any later version.

//

// This program is distributed in the hope that it will be useful,

// but WITHOUT ANY WARRANTY; without even the implied warranty of

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

// GNU General Public License for more details.

// https://www.gnu.org/licenses/



#pragma once



#include "defs.h"

#include "lbm.h"

#include "coefs.h"

#include "taylor.h"

#include "cs.h"



struct CORE : DEF

{

	LBM lbm;



	int ORDER=2;	// Taylor order

	int collision_model=0;

	int feqorder=-1; // negative = symbolic

	int verbosity=0;

	bool constvel=false;

	bool bigprint;

	bool EFDEprint;

	bool use_tay_vels=true;

	bool use_tay_meqs=true;

	bool use_tay_ms=true;

	bool use_symbolic_B=false;

	bool factor_coefs=true;



	char target_dir[1024];



	matrix I, E[MQ], A, B, C, M, Mi, BMi, e[MQ];

	matrix *As, *dBMi;

	vector<CS> vecCS;

	int CC[MDIM][MQ];

	matrix *CP;

	ex *trP;

	bool allocated=false;



	ex ***coefs_efde, ****coefs_epde;

	ex **sym_vels, **sym_fs, **sym_ms, **sym_feqs, **sym_meqs;

	// central values:

	ex *sym_m, *sym_meq, *sym_f, *sym_feq;



	TAYLOR **tay_vels, **tay_ms, **tay_meqs;



	matrix const_feq, const_sub_feq, const_meq;

	matrix *feq, *meq;

	bool Feq_allocated=false;



	ex sym_cs, sym_vel[MDIM];



	void taylor_list_vels(lst &list);

	void taylor_list_fs(lst &list);

	void taylor_list_ms(lst &list);

	void taylor_list_meqs(lst &list);





	void desc2str(char*str, int k, int j1, int j2, int j3);

	void number2str(char*str, const char*s1, const char*s2, int number);

	void coord2str(char*str, int x, int y, int z);

	void desc2str(char *s_space_coeff, char* s_time_coeff, char* s_space_formula, char* s_time_formula, int l, int i1, int i2, int i3);



	std::string ReplaceString(std::string subject, const std::string& search, const std::string& replace);

	inline void simplify_matrix(matrix &A) { simplify_matrix(A,Q,Q); }

	inline void simplify_matrix(matrix &A, int Nrows, int Ncols) { for (int i=0;i<Nrows;i++) for (int j=0;j<Ncols;j++) A(i,j) = expand(A(i,j)); }

	inline void empty_matrix(matrix &A) { empty_matrix(A,Q,Q); }

	inline void empty_matrix(matrix &A, int Nrows, int Ncols) { for (int i=0;i<Nrows;i++) for (int j=0;j<Ncols;j++) A(i,j)=0;; }



	inline int indO(int a, int b=0, int c=0) { return a + b*(ORDER+1) + c*(ORDER+1)*(ORDER+1); }

	void Meq_symbolic(matrix &M, int l, int x, int y, int z);

	void Feq_symbolic(matrix &F, int l, int x, int y, int z);

	void Feq_symbolic_Taylor(matrix &F, int l, int x, int y, int z);

	void fill_sym_vels_to_basis(COEFS &CF);



	void define_matrix_A_worker(matrix &_A, ex &u, ex &v, ex &w);

	void define_matrix_A();

	void define_matrix_B();

	void define_matrix_M();

	void define_equilibrium();

	void allocate();



	void compute_coefficients();

	void compute_EFDE_coefs();

	void compute_EPDE_coefs();



	void export_EFDE_coefs_bigprint(ex &c, ofstream &gout, char* str, bool last);

	void export_EFDE_coefs(int pde);

	void export_EPDE_coefs(int pde);

	void export_EPDE_coefs_bigprint(ex &expr, ofstream &fout, char *str);



	CORE(int iDIM, int iQ) : DEF(iDIM, iQ) { sprintf(target_dir,"default"); };

	~CORE();

};