test_solid_mechanics_model_reassign_material.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Jun 20, 09:48

test_solid_mechanics_model_reassign_material.cc
View Options

	/**
	* @file test_solid_mechanics_model_reassign_material.cc
	* @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
	* @date Mon May 6 17:37:17 2013
	*
	* @brief test the function reassign material
	*
	* @section LICENSE
	*
	* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
	* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
	*
	* Akantu is free software: you can redistribute it and/or modify it under the
	* terms of the GNU Lesser General Public License as published by the Free
	* Software Foundation, either version 3 of the License, or (at your option) any
	* later version.
	*
	* Akantu 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 Lesser General Public License for more
	* details.
	*
	* You should have received a copy of the GNU Lesser General Public License
	* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
	*
	*/

	/* -------------------------------------------------------------------------- */
	#include "static_communicator.hh"
	#include "solid_mechanics_model.hh"
	#include "material.hh"
	#include "aka_grid_dynamic.hh"
	using namespace akantu;


	class StraightInterfaceMaterialSelector : public MaterialSelector {
	public:
	StraightInterfaceMaterialSelector(SolidMechanicsModel & model,
	const std::string & mat_1_material,
	const std::string & mat_2_material,
	bool & horizontal,
	Real & pos_interface) :
	model(model),
	mat_1_material(mat_1_material),
	mat_2_material(mat_2_material),
	horizontal(horizontal),
	pos_interface(pos_interface) {
	Mesh & mesh = model.getMesh();
	UInt spatial_dimension = mesh.getSpatialDimension();

	/// store barycenters of all elements
	mesh.initElementTypeMapArray(barycenters, spatial_dimension, spatial_dimension);

	for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
	GhostType ghost_type = *gt;
	Element e;
	e.ghost_type = ghost_type;

	Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
	Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
	for(; it != last_type; ++it) {
	UInt nb_element = mesh.getNbElement(*it, ghost_type);
	e.type = *it;
	Array<Real> & barycenter = barycenters(*it, ghost_type);
	barycenter.resize(nb_element);

	Array<Real>::iterator< Vector<Real> > bary_it = barycenter.begin(spatial_dimension);
	for (UInt elem = 0; elem < nb_element; ++elem) {
	e.element = elem;
	mesh.getBarycenter(e, *bary_it);
	++bary_it;
	}
	}
	}
	}

	UInt operator()(const Element & elem) {
	UInt spatial_dimension = model.getSpatialDimension();
	const Vector<Real> & bary = barycenters(elem.type, elem.ghost_type).begin(spatial_dimension)[elem.element];

	/// check for a given element on which side of the material interface plane the bary center lies and assign corresponding material
	if (bary(!horizontal) < pos_interface) {
	return model.getMaterialIndex(mat_1_material);;
	}
	return model.getMaterialIndex(mat_2_material);;

	}

	bool isConditonVerified() {

	/// check if material has been (re)-assigned correctly
	Mesh & mesh = model.getMesh();
	UInt spatial_dimension = mesh.getSpatialDimension();
	GhostType ghost_type = _not_ghost;

	Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
	Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
	for(; it != last_type; ++it) {
	Array<UInt> & el_idx_by_mat = model.getElementIndexByMaterial(*it, ghost_type);
	UInt nb_element = mesh.getNbElement(*it, ghost_type);
	Array<Real>::iterator<Vector<Real> > bary = barycenters(*it, ghost_type).begin(spatial_dimension);
	for (UInt elem = 0; elem < nb_element; ++elem, ++bary) {
	/// compare element_index_by material to material index that should be assigned due to the geometry of the interface
	UInt mat_index;
	if ((*bary)(!horizontal) < pos_interface)
	mat_index = model.getMaterialIndex(mat_1_material);
	else
	mat_index = model.getMaterialIndex(mat_2_material);

	if (el_idx_by_mat(elem,0) != mat_index)
	/// wrong material index, make test fail
	return false;


	}
	}
	return true;
	}

	void moveInterface(Real & pos_new, bool & horizontal_new) {
	/// update position and orientation of material interface plane
	pos_interface = pos_new;
	horizontal = horizontal_new;
	model.reassignMaterial();
	}

	protected:
	SolidMechanicsModel & model;
	ElementTypeMapArray<Real> barycenters;
	std::string mat_1_material;
	std::string mat_2_material;
	bool horizontal;
	Real pos_interface;
	};


	/* -------------------------------------------------------------------------- */
	/* Main */
	/* -------------------------------------------------------------------------- */
	int main(int argc, char *argv[]) {

	bool test_passed;

	debug::setDebugLevel(dblWarning);
	initialize("two_materials.dat", argc, argv);

	/// specify position and orientation of material interface plane
	bool horizontal = true;
	Real pos_interface = 0.;

	UInt spatial_dimension = 3;

	akantu::StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator();
	akantu::Int psize = comm.getNbProc();
	akantu::Int prank = comm.whoAmI();

	Mesh mesh(spatial_dimension);

	akantu::MeshPartition * partition = NULL;
	if(prank == 0) {
	/// creation mesh
	mesh.read("cube_two_materials.msh");

	partition = new akantu::MeshPartitionScotch(mesh, spatial_dimension);
	partition->partitionate(psize);
	}

	/// model creation
	SolidMechanicsModel model(mesh);
	model.initParallel(partition);
	delete partition;

	/// assign the two different materials using the StraightInterfaceMaterialSelector
	StraightInterfaceMaterialSelector * mat_selector;
	mat_selector = new StraightInterfaceMaterialSelector(model, "mat_1", "mat_2", horizontal, pos_interface);

	model.setMaterialSelector(*mat_selector);
	model.initFull(SolidMechanicsModelOptions(_static));
	MeshUtils::buildFacets(mesh);
	// model.setBaseName("test_reassign_material");
	// model.addDumpField("element_index_by_material");
	// model.addDumpField("partitions");
	// model.dump();

	/// check if different materials have been assigned correctly
	test_passed = mat_selector->isConditonVerified();
	if (!test_passed) {
	AKANTU_DEBUG_ERROR("materials not correctly assigned");
	return EXIT_FAILURE;
	}

	/// change orientation of material interface plane
	horizontal = false;
	mat_selector->moveInterface(pos_interface, horizontal);

	// model.dump();

	/// test if material has been reassigned correctly
	test_passed = mat_selector->isConditonVerified();
	if (!test_passed) {
	AKANTU_DEBUG_ERROR("materials not correctly reassigned");
	return EXIT_FAILURE;
	}

	finalize();

	if(prank == 0)
	std::cout << "OK: test passed!" << std::endl;

	return EXIT_SUCCESS;
	}

	/* -------------------------------------------------------------------------- */

test_solid_mechanics_model_reassign_material.ccNo OneTemporaryActions

File Metadata

test_solid_mechanics_model_reassign_material.ccView Options

Event Timeline

test_solid_mechanics_model_reassign_material.cc
No OneTemporary
Actions

test_solid_mechanics_model_reassign_material.cc
View Options